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SECTION  I 
INTRODUCTION 


NOVA-2' (Nuclear  Overpressure  Vulnerebility  Analysis,  Version  2)  is 
e digital  computer  program  representing  a complex  analysis  of  aircraft 
structural  elements  subjected  to  nuclear  overpressure  effects.  This 
version  contains  extensive  modifications  of  NOVA  (ref.  1,  now  also 
referred  to  as  NOVA-1)  with  expanded  analytical  capabilities,  improved 
accuracy,  and  more  efficient  computer  utilization.  This  document 
presents  a complete  description  of  the  analytical  methods  used  and  of 
the  computer  program.  Including  guidelines  for  running  the  program  and 
example  problems. 

The  NOVA-2  program  provides  a Technique  for  predicting  the  elastic 
and  inelastic  response  of  aircraft  structural  elements  to  the  transient 
pressure  loads  associated  with  the  blast  wave  from  a nuclear  explosion. 
These  high  intensity  pressure  loads  are  treated  separately  from  the  gust 
loads  due  to  the  blast  wave  and  are  associated  with  the  Initial  reflec- 
ted pressure  which  occurs  curing  diffraction  of  the  blast  wave  around 
the  structure.  Because  the  pressures  exist  for  such  a short  time,  they 
excite  high  frequency,  secondary  structure  such  as  skin  panels,  string- 
ers, longerons,  frames,  ribs,  canopies  and  radomes.  The  gust,  or  post 
diffraction,  loads  tend  to  excite  low  frequency,  primary  structural 
surfaces  such  as  the  wings,  fuselage,  and  horizontal  and  vertical  tails. 
Therefore,  the  separation  between  pressu  *e  (often  referred  to  as  over- 
pressure) and  gust  effects  on  aircraft  structure  is  generally  based  on 
secondary  high  frequency  structure  and  primary  low  frequency  structure. 

The  program  determines  the  slant  range  between  the  aircraft  and 
point  of  burst  for  specified  levels  of  structural  damage  as  a function 
of  structural  element;  weapon  yield,  orientation  and  height  of  burst; 
aircraft  speed  and  altitude;  and  degree  of  probability  that  the  level  of 
damage  is  incurred  or  exceeded.  Likewise,  it  determines  the  response  of 
structural  elements  for  a specified  slant  range. 
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A single  element  dynamic  analysis  technique,  which  considers  both 
linear  elastic  and  inelastic  deformations  and  assumes  that  the  element 
does  not  interact  with  adjacent  elements,  reduces  the  complexity  of  the 
modeling  and  analysis,  and  thus  provides  a solution  more  rapidly  than  a 
finite  element  analysis.  The  errors  introduced  by  considering  isolated 
structural  elements  rather  than  multi  elements  appear  to  be  reasonable 
for  most  components  when  compared  with  those  introduced  by  other  uncer- 
tainties associated  with  the  analysis.  However,  for  structures  with 
rapidly  changing  cross  section,  such  as  a nose  radome,  the  errors  may  be 
unacceptable. 

The  program  uses  a quasi-strip  method  for  predicting  the  pressure 
loading  on  the  lifting  surfaces  of  the  aircraft  prior  to  the  arrival  of 
the  blast  wave.  The  pressure  induced  by  the  thickness  of  the  airfoil  is 
neglected  and  the  pressure  loading  is  assumed  to  be  a function  of  and 
directly  proportional  to  the  angle  of  attack.  The  pressure  loadings 
induced  on  the  wing  and  tail  surfaces  by  the  blast  wave  are  predicted 
using  linear  acoustic  theory  and  the  assumption  that  the  airfoil  sec- 
tions can  be  represented  by  their  mean-camber  lines. 

The  program  predicts,  prior  to  and  during  diffraction  of  the  blast 
wave,  the  pressure  loading  on  the  surface  of  the  fuselage  by  relating 
the  fuselage  to  an  equivalent  body  of  revolution  with  the  same  dis- 
tribution of  cross-sectional  area  along  its  length  and  by  applying  the 
reflection  theory  for  the  interaction  of  a shock  wave  with  a flat 
surface. 

A 1-KT  nuclear  standard,  based  on  data  obtained  from  the  AFWL 
SPUTTER  and  SAP  fluid  dynamics  programs,  provides  the  time-dependent 
free-air  blast  characteristics  for  the  BLAST  routines.  For  near-ground 
bursts  where  the  blast  wave  strikes  the  ground  and  is  reflected,  two 
models  are  now  optional  in  NOVA-2:  1)  an  empirical  model,  which  describes 

the  wave  forms  associated  with  regular  reflection  and  the  transition  to 
Mach  reflection,  and  2)  a comprehensive  data  tape,  based  on  the  REFLECT 
code,  which  provides  a time  history  of  the  ground-reflected  blast  char- 
acteristics. 


The  program  consists  of  three  distinct  routines,  NOVA,  DEPROB 
(Dynamic  Elastic  Plastic  Response  of  Beams),  and  DEPROP  (Dynamic  Elastic 
Plastic  Response  of  Panels),  written  in  FORTRAN  IV. 

The  NOVA  routine  is  the  master  routine  which  controls  the  logic  of 
the  overall  program.  It  contains  the  subroutines  for  predicting  the 
aerodynamic  flight  loads  and  the  blast  pressure  loads  that  are  applied 
to  the  lifting  surfaces  and  fuselage  during  subsonic  and  supersonic 
flight,  and  for  determining  the  slant  range  at  wh£ch  a structural  ele- 
ment Incurs  damage  which  has  been  specified  on  a probabilistic  basis. 

The  DEPROB  routine  provides  the  response  of  aircraft  structure  such 
as  stringers,  longerons,  frames,  ribs,  and  conical  or  cylindrical 
radomes  which  can  be  represented  by  an  annular  cross  section.  The 
method  of  analysis  used  in  this  routine  applies  to  beams  which  can  be 
modeled  in  one  dimension  by  a series  of  discrete  masses  interconnected 
by  weightless  bars.  Major  additions  to  DEPROB  in  NOVA-2  are  the  ability 
to  analyze  elements  with  variable  cross  section,  the  addition  of  simply 
supported  and  free  edge  conditions,  an  Improved  elastic-plastic  stress- 
strain  model,  and  the  inclusion  of  rib  buckling  as  a failure  mechanism. 

The  DEPROP  routine  provides  the  response  of  aircraft  skin  panels, 
canopies,  and  radomes  that  can  be  approximated  by  a cylindrical  panel. 
The  linear  elastic  option  applies  to  single  and  multilayered  panels  of 
isotropic  or  orthotropic  material,  and  the  elastic-plastic  option 
applies  to  single  layer  panels  of  Isotropic  material.  DEPROP  has  been 
modified  to  include  the  following:  symmetric  or  nonsymmetric  combina- 

tions of  clamped  or  simply  supported  edge  constraints,  a much  improved 
eltstic-plastic  stress-strain  model,  and  Improved  overall  accuracy. 

The  NOVA-2  program  as  a whole  represents  a much  more  efficient 
code  than  NOVA-1,  requiring  less  computer  core  and  time.  This  is 
accomplished  by  making  more  extensive  use  of  overlaying  techniques  and 
through  improved  computational  techniques. 
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SECTION  II 
BLAST  MODELS 


In  order  co  simulate  the  overpressure  effects  of  a blast  wave 
emanating  from  a nuclear  explosion,  NOVA-2  provides  the  analyst  with  a 
choice  of  ground  reflection  models  when  ground  reflection  is  important. 
A one-dimensional  free-air  blast  model  is  the  basis  of  both  models,  and 
is  used  when  ground  reflection  is  not  important.  One  ground  reflection 
model  is  a semiemplrlcal  functional  curve  fit  (that  used  in  NOVA-1) ; 
the  second  is  a two-dimensional  data  tape  generated  by  a hydrodynamic 
computer  code.  These  two  models  are  described  in  subsections  2.1  and 
2.2,  respectively. 

The  free-air  blast  model  in  NOVA-2  is  the  most  recent  AFWL  1-KT 
curve  fit  model  available  (refs.  2 and  3),  which  is  based  on  the  AFWL 
1-KT  nuclear  standard  data  tape  (ref.  4).  Figure  1 compares  overpres- 
sure versus  range  for  both  free-air  models.  For  historical  complete- 
ness, the  IBM  Problem-M  curve  contained  in  reference  5 is  also  shown. 

By  coincidence  this  curve  nearly  coincides  with  that  of  the  current 
curve  fit  model,  and  the  differences  arc  indistinguishable  in  figure  1. 


2.1  ANALYTICAL  GROUND-REFLECTED  BLAST  MODEL 

The  first  ground-reflected  blast  model  is  a semiempirical,  two- 
dimensional  model,  identical  in  form  to  that  used  in  NOVA-1.  The  only 
differences  are  minor  changes  in  curve  fit  coefficients.  The  basic 
mathematical  model  is  documented  in  reference  6. 


2.2  REFRA  GROUND-REFLECTED  BLAST  MODEL 

The  second  ground  reflected  blast  model  available  in  NOVA-2  is  in 
the  form  of  a two-dimensional  blast  tape  (logical  file  TAPE10) , read  and 
interpreted  by  the  REFRA  routine  (ref.  7).  The  original  data  base  is 
generated  by  the  REFLECT  code  (ref.  8)  and  subsequently  processed  for 
more  efficient  use  by  REFRA. 
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Only  data  pertaining  to  region  III  o£  figure  2 are  provided  on 
tape,  aa  the  undisturbed  region,  region  I,  merely  represents  ambient 
conditions,  and  data  for  region  II  are  provided  by  the  free-air  model. 

A moving  mesh  model  in  the  REFLECT  code  provides  data  behind  the  shock 
fronts  which  are  precisely  defined  and  free  of  the  complications  asso- 
ciated with  various  artificial  smoothing  processes.  The  REFRA  routine 
searches  the  data  tape  to  determine  which  region  is  appropriate  and 
returns  the  corresponding  blast  data.  It  should  be  readily  apparent 
that  because  the  data  tape  is  limited  to  regions  behind  the  Mach  and 
reflected  shocks,  great  savings  in  tape  length  and  computer  search  time 
are  realized.  Accordingly,  a very  comprehensive  time  history  of  the 
reflected  wave,  and  the  region  behind  it,  is  made  available  on  tape. 

Each  REFLECT  run,  and  hence  each  corresponding  data  tape,  is  char- 
acterized by  a unique  scaled  height  of  burst  (above  ground  level) , 
scaled  to  1-KT  at  sea  level.  To  select  the  appropriate  tape  for  the 
problem  of  interest,  figure  3 shows  the  relationship  of  relevant  geo- 
metrical parameters. 

Once  the  ground  altitude  (H  ),  the  aircraft  altitude  (H)  and  the 

S 

vertical  separation  between  the  aircraft  and  the  burst  (z)  are  selected, 

the  height  of  burst  is  uniquely  specified  (HOB  - H - H - z).  The 

8 

scaled  height  of  burst  (SHOB)  is  then  determined  from  the  following 
relationship,  making  use  of  modified  Sachs  scaling: 


SHOB  - 


w 


HOB 


where  W is  the  yield  in  kilotons,  p is  the  ambient  pressure  at  the  air- 
craft  altitude,  and  pQ  is  the  ambient  pressure  at  sea  level,  taken  to  be 
14.696  psi. 

It  can  be  seen  that  the  altitude  of  the  aircraft  relative  to  the 
burst  must  remain  constant  for  problems  using  a single  data  tape. 

Hence,  when  the  program  iterates  to  find  a critical  slant  range,  the 
iteration  is  restricted  to  constant  aircraft  and  burst  altitudes  for 
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Figure  3.  Geometrical  Relationship  of  Aircraft,  3urst  Center, 
and  Ground 
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this  ground  reflection  model,  whereas  the  method  described  in  subsec- 
tion 2.1  does  not  have  this  limitation  since  the  empirical  model  will 
handle  any  height  of  burst.  The  REFRA  model  will,  in  general,  also  take 
slightly  more  computer  time  for  problems  Involving  ground  reflection. 

It  should  be  noted,  however,  that  the  REFRA  model  offers  a much  more 
sophisticated,  complete  description  of  the  phenomena  of  ground  reflec- 
tion. 


SECTION  III 

PREBLAST  AND  BLAST- INDUCED  LOADING 


3.1  LOADING  ON  LIFTING  SURFACES 

Tha  pressure  loading  at  a point  on  a lifting  surface  of  an  air- 
craft engulfed  by  a blast  field  may  be  approximated  In  terms  of  two 
pressure  distributions;  the  first  associated  with  the  flow  prior  to 
blast  arrival  and  the  second  induced  by  the  blast  wave.  The  form  used 
is  as  follows: 


(AC  w 
p ss 


+ Ac  w ) t pV 
P g 2 r 


(1) 


where 

p.  and  p are  the  pressures  at  a point  on  the  lower  and  upper 

x»  u 

surfaces,  respectively 

ACp  and  Acp  are  the  steady-state  and  transient  pressure  coeffi- 
cients, respectively,  at  the  point  on  the  lifting 
surface 

w and  w are  the  steady-state  and  blast-induced  flow  velocities 
ss  g ' 

normal  to  the  lifting  surface,  respectively,  at  the 
point  on  the  surface 

p is  the  instantaneous  density  at  the  point  on  the 

surface 

is  the  resultant  velocity  of  the  steady-state  flow 
plus  the  blast-induced  flow  at  the  point  on  the 
surface 


The  procedures  by  which  AC^  and  Ac^  are  calculated  for 
in  the  analysis  are  presented  in  the  remainder  of  this 


lifting  surfaces 
subsection. 
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The  pressure  increase  on  the  lower  surface,  p,  - pa,  is  taken 
equal  in  magnitude  and  opposite  In  sign  to  the  pressure  decrease  on 
the  upper  surface,  p , - pm,  where  p^  Is  the  pressure  at  the  point  in 
the  absence  of  the  lifting  surface. 

h m p.  + \ (Pt  - pu) 

Pu  * P-  - 7 <Pg  - Pu)  <3) 

Since  p.  and  p cannot  physically  achieve  negative  values,  these  quan- 
x*  u 

titles  are  set  equal  to  zero  within  the  computer  analysis  whenever 
either  becomes  negative  according  to  equation  (2)  or  equation  (3) ; how- 
ever, the  correct  value  of  p.  - p , as  given  by  equation  (1),  is 

x>  u 

retained  In  the  analysis. 

3.1.1  Preblast  Loading 

This  section  provides  formulations  for  predicting  the  pressure 
loading  coefficient  ACp  on  lifting  surfaces  before  exposure  to  blast 
environment.  Although  more  elaborate  techniques  are  available,  l.e., 
surface  methods,  the  emphasis  here  is  on  practical  and  simplified  tech- 
niques which  are  more  appropriate  for  a vulnerability  code.  The  form- 
ulations presented  here  are  for  subsonic,  transonic,  and  "low  to 
medium"  supersonic  flight  ranges  and  consider  typical  classes  of  wing 
or  tall  planforms  during  a symmetric  flight  maneuver.  The  derived 
equations  form  the  basis  for  the  program  in  the  vulnerability  code. 

A rigorous  approach  cannot  be  offered  for  calculating  the 
pressure  coefficient  at  arbitrary  points  on  all  possible  planforms 
over  a broad  range  of  Mach  numbers.  From  the  practical  point,  it  is 
necessary  to  introduce  many  simplifications.  These  are  covered  in 
the  discussions  to  follow. 
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The  surfaces  ere  assumed  to  be  flat  plates,  This  means  that 
the  incremental  pressure  at  any  point  on  the  planfonn  is  due  to  angle  of 
attack  only,  and  thickness-induced  pressures  are  neglected.  This  asauop 
tion  is  a reasonable  one,  except  possibly  at  points  near  the  leading 
edge.  Furthermore,  the  angle  of  attack,  a,  is  assumed  to  be  small,  so 
that  the  pressure  distributions  are  directly  proportional  to  the  angle 
of  attack. 

The  general  planform  to  be  considered  is  depicted  in  figure  4. 
The  planform  is  assumed  to  have  a line  of  symmetry  which  is  parallel  to 
the  free-stream  direction.  The  tip  chord  is  approximated  by  a line 
parallel  to  the  centerline  of  the  planform. 


s 


8.  + S_ 

L T 


(4) 


For  purposes  of  calculating  the  pressure  coefficient,  ACp,  the  actual 
planform  is  replaced  by  a swept  wing  which  has  at  station  y the  same 


coordinates  and  slopes  of  the  leading  and  trailing  edgesj 

rdxi 


V |_dn 
figure  5 


^ — 


r r^Li 

V [dn  J ’ 


as  the  actual  planform.  Consider  the  specific  planform  in 
The  actual  (semispan)  planform  is  that  described  by  the 
straight  line  segments  OabcdeO.  In  calculating  the  pressure  at  point  1, 
the  equivalent  planform  will  be  OagfO.  On  the  other  hand,  for  point  2, 
the  equivalent  swept  wing  will  be  idchi.  Thus,  the  equivalent  swept 
wing  is  one  with  leading  and  trailing  edge  sweeps  equal  to  the  leading 
and  trailing  edge  sweeps  of  the  actual  wing  at  station  y.  This  method 
may  be  considered  as  a "quasi-strip"  method  in  that  the  strip  at  the 
pressure  spanwise  station  determines  the  planform.  Yet  it  is  not  a 
strip  method  per  se  because  finite  span  effects  are  accounted  for 
approximately.  In  fact,  if  the  slopes  of  the  leading  and  trailing  edges 
of  the  actual  planfonn  are  constant,  the  equivalent  planform  is  the 
actual  planfonn  regardless  of  where  the  pressure  is  to  be  determined. 
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Figure  4.  Actual.  Planform  Geometry 
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Two  possibilities  for  the  semispan,  s,  can  arise:  Referring 
to  figure  6,  the  semi-span  of  the  equivalent  wing  can  be  s with  the  tip 
chord  greater  than  zero  or  s^  with  no  tip  chord.  The  situation  in  (b) 
will  not  occur  too  often  for  actual  wings.  An  example  when  this  occurs 
is  given  by  (c),  where  the  pressure  is  to  be  computed  for  a point 
y < y.  Provisions  are  made  in  the  formulation  and  in  the  computer 
program  to  handle  cases  of  types  (b)  and  (c),  although  they  are  not 
expected  to  be  used  very  often.  These  cases  are  numbered  2,  7,  8 and 
9 among  the  possibilities  listed  in  table  1. 

The  equivalent  planform  geometry  is  presented  in  figure  6. 

The  coordinate  systems  (x,  y)  or  (€>n)  of  the  original  planform  are 
retained  for  convenience,  so  that  the  coordinate  of  the  pressure  point 
(x,y)  on  the  actual  planform  is  also  the  (x,y)  on  the  equivalent  plan- 
form.  It  should  be  noted  that  the  leading  and  trailing  edges  Intersect 
the  centerline  at 


*1,(0) 

* *l  - y tan^ 

(5) 

xT(0) 

- xT  - y tanAT 

(6) 

When  the  leading  or  trailing  edge  is  curved  or  is  composed  of 
segments  of  straight  lines,  several  equivalent  planforms  are  possible. 
Each  such  planform  will  have  its  leading  and  trailing  edges  described 
by 


x^n)  * xL  + (n  - y)  tanAL  (7) 

x^n)  * xT  + (n  - y)  tanAT  (8) 

and,  in  general,  the  planform  apex  will  not  be  at  (£,n)  ■ (0,0).  The 
chord  at  n is  defined  by  the  expression 


-15- 


r 


PUS  I mil  II  UW,  I " I -T-H  I w">  ’I  — 


c(n)  - x^n)  - x^n)  ■ (xT  - x^  + (n-y) (tanAT-tanAL)  (9) 


The  chord  becomes  zero  when 


(10) 

If  s^  > a,  Che  equivelenc  plenform  ie  terminated  at  n ■ 8,  and  corre- 
sponds to  (a),  figure  6. 

If  s^  < s,  the  equivalent  plenform  is  terminated  at  n ■ s^  and  the 
taper  ratio  will  be  zero.  This  corresponds  to  (b)(  figure  6. 


y + 


*T  - *L 


ta 


- 


tanA. 


= i, 


Referring  to  table  1,  Cases  1 and  2 are  subsonic  situations, 
Cases  3-9  are  supersonic,  while  Case  10  is  transonic.  For  supersonic 
flight,  distinctions  have  to  be  made  to  reflect  leading  and  trailing 
edge  conditions.  Let  B - If  6 cot  A > 1,  the  edge  is  super- 

sonic; if  6 cot  A < 1,  the  edge  is  subsonic. 


In  the  absence  of  a better  way  of  determining  transonic  pres- 
sures, Case  10  is  programmed  to  proceed  as  follows:  If  0.8  < M < 1.2, 

the  pressure  is  computed  for  M ■ 0.8  and  M ■ 1.2  according  to  appro- 
priate Cases  1-9.  The  pressure  is  then  interpolated  linearly  for  the 
pressure  coefficient  at  the  desired  Mach  number  according  to 


j 

1 

i 


(ACp)0.8<M<1.2 


(Hjr)  <4CP Ws + (Srr4) 


(ID 


Table  2 presents  a list  of  certain  planform  parameters  which 
are  useful  in  the  development  of  the  pressure  coefficients. 


Different  approaches  are  followed  for  the  subsonic  and  super- 
sonic  cases.  For  supersonic  cases,  the  pressures  may  be  obtained 
directly  from  linearized  steady  supersonic  flow,  with  additional 


-18- 


f 


PERTINENT  PLANFORM  PARAMETERS 


of  reference  9 is  assumed  to  be  equal  to  1.0  and 


approximations  for  certain  regions  of  the  planform.  In  contrast,  for 
the  subsonic  cases,  the  spanwlse  loadings  (lift  per  unit  span)  are 
first  determined.  The  pressures  are  then  determined  by  assuming  a 
chordwise  pressure  distribution  and  determining  an  "amplitude"  for  this 
distribution  such  that  the  pressures  integrated  over  the  chord  give 
the  correct  lift  per  unit  span  at  that  station.  The  distribution  used 
for  all  subsonic  Mach  numbers  is  a modification  of  the  chordwise  dis- 
tribution from  thin  airfoil  theory,  and  is  illustrated  in  figure  7. 

Each  of  the  cases  will  now  be  discussed  individually. 


a.  Subsonic  Cases  (M  < 0.8) 

Given  the  location  (x,y)  on  the  planform,  the  quantities 
x^,  xT,  tanA^,  tanA^  can  be  determined  from  the  planform  geometry.  If 

*T  ' \ 


8 ■ v 4* 

1 1 tanAr  - tanA 


> s.  Case  1 is  considered;  if  < s.  Case  2 


L T 

is  appropriate.  Only  Case  1 will  be  discussed  because  Case  2 is  iden- 
tical with  Case  1,  except  s is  replaced  by  s^  and  the  taper  ratio  A 
is  set  equal  to  zero.  The  parameters 


i ■ *’  V x- BA  (-lr) 

may  be  computed  according  to  the  expressions  given  in  table  2.  Using 
these  parameters,  the  spanwlse  loading  at  y 


(12) 


k of  reference  9 


is  taken  equal  to  1.0. 
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and 


ecT 


ecT 


, since  < - 1 


(13) 


may  be  obtained  by  interpolation  from  the  prepared  curves  of  reference  9. 
Since  these  interpolations  are  over  multiple  variables,  expressions 
were  derived  for  approximating  the  curves  of  figures  3 and  4 of  ref- 
erence 9 to  facilitate  the  interpolation  process.  The  pressure  differ- 
ence coefficients  are  then  computed  according  to 


where 


1CP  ■ \ * i (i)  f<*>  (2 


+ 

XT  - *L 


[tan  A - tan 


v) 


(14) 


11 


-Pi)- 


xL  x < 0.05  xT  + 0.95  Xj^ 


f(x) 


(15) 


.6.1287 


lX T ' X 
V X - x. 


- , 0.05xT+0.95xL  £ x < xT 


b.  Supersonic  Cases  (M  > 1.2) 

For  supersonic  flow,  the  planform  is  divided  into  distinct 
regions  by  Mach  lines  emanating  from  leading  edges  of  the  root  chord  and 
the  tip  chord  and  by  "reflected"  Mach  lines  from  other  edges  as  shown 
in  figure  8.  For  each  of  these  regions,  different  pressure  expressions 
apply.  The  particular  region  in  which  the  point  of  interest  (x,y)  lies 
must  be  identified  as  to  type  of  region  so  that  an  appropriate  pressure 
formula  may  be  used.  To  illustrate  this  brief  description,  consider 
Case  3.  Depending  on  the  sweep  of  the  Mach  line,  A^  * tan  ^S,  three 
possible  situations  may  arise.  These  situations  are  depicted  in 
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figure  8 and  are  described  in  references  10  and  11.  For  points  in 
Region  1,  the  AC^  would  be  the  same  ca  on  an  infinite  wing  with  the 
same  sweep  as  the  leading  edge,  and  is 


4 cot  A. 


(AC  ) 


P'l 


{ 


(16) 


2 2 

g cot  Al  - 1 


For  Region  II,  the  AC  may  be  obtained  from  the  results  for  the  corre- 
sponding region  on  a delta  wing  with  supersonic  leading  e-’ges. 


4 cot  A. 


(AC) 


PH 


V 


2 2 

8 cot  A.  - 1 
L 


1 - ^ sin-1  1 1 - >2  tan2  v 

n '22  2 

8 (cot  A^-tan  v) 


i 

i 


(17) 


where 


„ 2 
tan  v 


. I 


| x ' (\  ' yL  tan  AL}  j 


For  Region  III,  the  AC^  may  be  defined  by: 


<Vm 


8 cot  A 


tan 


-1 


J (1  + 

i cot  Al[x-(xl- 


3 cot  A^) (s-y) 


y tan  AL)-s(tan  A^)-6(s-y)] 


(18) 
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Simple  expressions  cannot  be  found  for  Regions  IV  and  V.  For  a point 

(x,y)  in  Region  IV  of  figure  8(b),  the  ACp  is  interpolated  linearly 

between  the  pressures  at  points  a and  b where  (AC^)^  and  (ACp)^  apply. 

If  (x,y)^  lies  in  IV,  but  ahead  of  the  dashed  line  in  figure  8(c),  the 

same  procedure  may  be  applied.  For  (x,y)^  in  IV  but  behind  the  dashed 

line,  the  ACp  is  interpolated  linearly  between  points  c and  d;  point  c 

has  (ACp)  and  ACp  for  point  d is  exactly  0.  For  points  in  Region  V, 

such  as  (x,y)_,  AC  is  also  interpolated  linearly  between  points  c and 
2 P 

d.  In  all  cases,  the  line  along  which  interpolation  takes  place  is 
parallel  to  the  trailing  edge. 

The  Mach  line  patterns  for  Cases  4,  5,  and  6 are  as  shown  in 
figure  9.  For  the  above  cases,  only  (AC^)y^  can  be  obtained  from  known 
solutions.  Using  the  solution  for  a delta  wing  with  subsonic  edges 
(refs.  10  and  11),  (ACp)Vj  is  defined  in  terms  of  the  present  coordinate 
system  by 


4 cot 


(AC  ) 


p'VI 


h 


2 2 

tan  Al  tan  v E(k) 


(19) 


where 


_ 2 
tan  v 


' \ * 


y tan 


VI 


and 


E(k)  ■ complete  elliptical  integral 


»/2  


2 2 

k sin  z dz;  k 


* 


2 2 
8 cot  At 


(20) 


For  Regions  VII  and  VIII,  etc.,  approximate  pressures  may  be  obtained 
through  Interpolation,  using  the  pressures  elong  the  boundaries  of 
Region  VI,  and  the  fact  that  the  pressure  must  go  to  zero  at  the  sub- 
sonic trailing  edge  and  the  tip  edge. 

Cases  7,  8 and  9,  which  are  shown  in  figure  10,  can  be 
considered  as  subcases  of  the  previous  Case  3,  and  exhibit  the  seme  type 
of  regions.  Therefore,  the  AC^'s  obtained  previously  also  apply  here. 

c.  Transonic  Cases  (0.8  < M < 1.2) 

The  pressure  coefficient  (ACp)  for  0.8  < M < 1.2  msy  be 
estimated  by  linear  interpolation  over  M utilizing  the  AC^'s  for 
M ■ 0.8  and  M ■ 1.2,  as  shown  by  equation  (11). 

3.1.2  Blast  Airloads  on  Lifting  Surfaces 
a.  General  Discussion  of  Method 

The  equations  for  the  aerodynamic  loads  coefficient, 

ACp  due  to  a blast  wave  are  presented  in  this  section.  These  coeffi- 
cients are  used  in  equation  (1)  for  determining  the  local  pressures. 

The  problem  to  be  addressed  is  the  determination  of  Ac ^ as  a function 
of  position  and  time  for  arbitrary  aircraft  speed,  blast  orientation, 
and  blast  strength. 

Equations  for  predicting  the  blast-induced  airloads  on 
lifting  surfaces  for  arbitrary  strength  of  the  blast  shock  are  not 
available.  References  12  through  18  have  however  demonstrated  that 
predictions  of  the  difference  in  loading  between  opposite  surfaces  of  a 
wing  or  tail  surface  exhibit  extensive  areas  of  good  agreement  with 
measured  loadings  where  the  prediction  is  based  on  acoustic  theory. 

Acoustic  theory  as  applied  to  thin  airfoils  is  based  on 
the  assumption  that  the  airfoil  section  of  a wing  or  tail,  for  example. 


11  ^Vl 


Is  chin  enough  ralativa  to  tba  chord  and  apan  that  tha  sactlon  can  ba 
raplacad  by  a line,  In  thla  cua  tha  caabar  lina.  Tha  flow  is  assuoad 
to  ba  attachad  to  tha  aurfaca.  Acoustic  thaory  has  been  demonstrated  by 
rafarencas  12  through  IS  to  provlda  fairly  good  agreement  with  measured 
airloads  due  to  blast  and  shock  wavas  whan  tha  comparison  vas  Bade  in 

tens  of  Ac  • There  are  a few  exceptions  to  this  which  oust  be  discussed. 

P 

One  particular  region  of  tha  shock  loading  that  has  been 

found  to  differ  froa  acoustic  predictions  is  in  tha  vicinity  of  tha 

shock  front  of  a blast- type  wave.  Tha  osaaurad  distribution  of  Ac  is 

P 

found  to  have  a flat-top  peak  hare,  in  contrast  to  tha  singularity 
predicted  by  acoustic  theory,  a.g.,  reference  17,  figure  5.  However,  it 
is  important  to  note  that  tha  integral  over  tha  airfoil  of  the  neasured 
loading  coefficient  Acp  is  found  to  agree  quite  wall  with  the  theoreti- 
cal prediction  of  acoustic  thaory,  even  to  blast-induced  angles  of 
attack  as  high  as  30  degrees,  e.g.,  reference  17,  figure  6a.  This 
result  is  similar  to  tha  well-known  leading-edge  singularity  for  sub- 
sonic airfoils  predicted  by  linearized  thaory,  where  in  practice,  Acp  is 
found  to  rise  to  large  values  near  tha  leading  edge.  Higher  observed 
values  of  Acp  further  rearward  of  tha  leading  edga  are  found  to  compen- 
sate for  tha  absence  of  any  singularity , however,  tending  to  make  the 
integral  over  tha  airfoil  of  neasured  Acp  agree  better  with  theoretical 
predictions. 

At  large  blast-induced  angles  of  attack,  the  flow  even- 
tually separates  from  the  upper  surface,  causing  Acp  to  drop  below  the 
predictions  of  thin-alrfoil  theory.  Experiments  which  demonstrate  the 
developstent  of  this  separation  with  time  are  described  in  references  12 
and  17.  Because  the  separation  affects  the  loading  at  long  times  rela- 
tive to  the  time  expected  for  a typical  aircraft  surface  to  respond 
structurally  to  the  blast  overpressure,  no  attempt  is  made  here  to 
include  the  complicated  variations  in  loading  history  at  the  late  times 
associated  with  flow  separation.  It  should  also  be  noted  that  the  sepa- 
ration occurs  on  the  low-pressure  side  of  a wing  or  tail  which  would 
generally  not  be  critical  to  the  overpressure  effects. 
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Ic  should  also  bs  noted  that  the  nonlinearities  in  the 


pressure  on  the  two  surfaces  tend  to  compensate  each  other;  the  nonlinear 
increase  in  - pa  with  angle  of  attack  is  accompanied  by  a nonlinear 
increase  in  p^  - p^  , although  the  latter  tends  to  be  smaller.  There- 
fore, ACp  tends  to  show  less  nonlinearity  with  angle  of  attack  than  do 
either  of  the  other  quantities.  No  general  technique  is  available  at 
present  for  predicting  nonlinear  loadings  for  the  general  cases  of 
Interest,  including  effects  of  shock  orientation,  airfoil  speed,  etc. 

The  linear  acoustic  theory  has  been  developed  for  a wide  variety  of 
cases,  so  it  is  employed  in  the  present  work. 

After  a blast  wave  impinges  on  a rectangular  wing  or 
tail,  the  pressures  generally  tend  to  return  to  equilibrium  much  more 
rapidly  from  waves  moving  in  a chordwise  direction  than  moving  spanwise. 
Therefore,  the  analysis  will  be  based  on  strip  theory,  considering 
strips  in  a chordwise  direction.  Spanwise  effects  are  expected  to 
take  place  over  times  which  are  long  compared  to  the  structural  response 
of  interest  for  overpressure. 

Table  3 contains  seven  blast  loading  cases  of  interest 

as  functions  of  the  airfoil  Mach  number,  M,  and  the  gust  Mach  number, 

M , of  the  blast  shock.  For  these  cases,  there  are  seven  forms  of  the 

g 

equation  for  Ac^. 

In  the  application  of  the  equations,  both  Mach  numbers 
are  taken  relative  to  the  flow  behind  the  blast  shock,  which  was  demon- 
strated in  reference  14  to  provide  the  best  correlation  of  the  airloads 
with  experimental  data.  The  blast  radius  is  assumed  to  be  very  large 
relative  to  the  dimensions  of  a wing  or  tail  chord,  so  the  blast  prop- 
erties could  be  taken  at  any  convenient  point  relative  to  the  chord; 
the  blast  properties  will  be  taken  at  the  first  point  of  blast  inter- 
cept, either  the  leading  or  the  trailing  edge  of  the  airfoil,  for 
determining  the  blast  Ac  . 
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BLAST  LOAOIMO  CASES  OF  INTEREST 


Th«  Mach  number  of  Che  airfoil  relative  to  the  flow 
behind  Che  shock  is  M - V^/a where  is  Che  velocity  of  the  fluid 
behind  Che  shock  relative  Co  Che  airfoil  and  a2  is  the  speed  of  sound 
behind  Che  shock.  The  velocity  diagram  is  shown  in  figure  11,  where  V 
is  the  airfoil  velocity  relative  to  the  undisturbed  fluid  ahead  of  the 
shock,  and  v^  is  the  fluid  velocity  of  the  blast  wave. 

Relative  to  the  fluid  behind  the  shock,  the  shock  velocity 
is  - v^,  where  V#  is  the  shock  velocity  relative  to  the  undisturbed 
fluid  ahead  of  the  shock. 

The  shock  correlations  carried  out  in  references  14  and 
17  indicate  that  the  envelopment  rate  of  the  shock  wave  is  to  be  con- 
sidered relative  to  the  "equivalent"  airfoil  sketched  in  figure  11, 
where  the  equivalent  airfoil  is  aligned  with  the  total  flow,  Vr,  behind 
the  blast  shock.  Then,  the  time  At  for  the  undisturbed  shock  to  pass 
over  the  airfoil  is  given  by 


V 

r 


v - v. 

s b 

cos($  - a2) 


(21) 


where  $ is  the  angle  between  the  shock  outward  normal  and  the  plane  of 
the  actual  airfoil.  The  positive  sign  applies  for  shocks  impinging  at 
the  leading  edge  first  and  the  minus  sign  for  shocks  impinging  at  the 
trailing  edge  first.  Defining  the  gust  Mach  number,  M^,  of  the  shock 
by 


M 

g 


cos  (*  “ 


(22) 


and  then  combining  the  equations  gives 


+ 


(23) 


M + M 

8 


a2 


c 

At 


Note  that  M can  be  poaitive  or  negative,  depending  upon  the  shock 

g 

orientation  relative  to  the  equivalent  airfoil;  positive  for  shocks 
arriving  from  the  front  and  negative  for  shocks  arriving  from  the  rear. 
Also,  because  a shock  wave  is  subsonic  relative  to  the  fluid  behind  it, 
the  absolute  value  of  M can  be  less  chan  unity,  thus  explaining  the 
presence  of  such  cases  in  table  3. 


b.  Discussion  of  Methods  for  Various  Cases  of  Airfoil  and 
Shock  Mach  Numbers 

The  wave  diagrams  used  here  for  the  seven  cases  are  drawn 
with  the  coordinates  fixed  on  the  undisturbed  fluid  at  rest  ahead  of  the 
blast  wave  in  the  absence  of  the  airfoil.  The  airload  equations  which 
correspond  to  figures  12  through  18  are  listed  in  tables  4 through  10. 
Acoustic  theory  for  thin  airfoils  assumes  that  all  perturbations  to  the 
fluid  due  to  the  airfoil  and  the  shock  wave  are  weak.  Therefore,  acous- 
tic theory  Is  based  on  the  assumption  that  the  angle  of  attack  at  all 
times  is  small  relative  to  a radian. 

The  X axis  is  directed  forward  with  the  origin  at  the 
point  of  shock  intercept  with  the  leading  edge  (or  with  the  trailing 
edge,  if  the  trailing  edge  is  intercepted  first)  of  the  airfoil.  The 
trace  of  the  leading  edge  is  represented  by  Line  I in  figures  12  through 
18.  The  X coordinate  is  scaled  with  the  chord,  c,  so  X ■ -1  at  the 
trailing  edge  at  the  time  of  intercept.  The  other  coordinate  in  the 
figures  is  reduced  time  T,  where  T is  equal  to  real  time  t scaled  with 
the  ambient  speed  of  sound,  a,  and  the  chord,  c.  The  trace  of  the 
trailing  edge  of  the  airfoil  is  represented  by  Line  II. 

The  equations  for  Ac^  are  given  in  terms  of  the  actual 
distance  x and  the  real  time,  with  x • c at  the  trailing  edge.  Actual 
distances  and  times  are  related  to  the  scaled  distances  and  time  by 
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Figure  16.  AirLoad  Regions  on  Airfoil  for 

Case  5 (M  > 1,  M > l) 

& 


-1  o 


Figure  L7.  AirLoad  Regions  on  AirfoiL  for 
Case  6 (M  > l,  0 *•  • 1) 


Table  4 


AIRLOAD  EQUATIONS  FOR  CASE  1 (M  < 1,  Mg  > 1) 


Line 

I 

II 

III 

IV 

V 

VI 

VII 


Equation 

0 
c 


(M  + Mg) at 
(M  + 1) at 
1 + M 


M + M„ 


c + (M  - l)at 


1 + M 


c + (M  - l)at 


Region  1 


2c 

( i+M) (1-M) a 


Ac  (t)  * — 

P1  [M*  - 1] 


|M„1 


M 


Region  2 


Ac  (x,  t)  » Ac 
P2  Pi 


1 - | tan*1 

' (M+Mg)  [ (l-M)at  +x]-(l+Mg)c  ‘ 

(1  + M ) (c  - x) 

Region  3 


Ac  ( x , t ) 
?3 


Ac 


1 - — tan"l 


M - Mg  (M  + l)at  - X 


M - 1 

l g 


1/2 


2 

T 


M 


Mg (1  + M) 


(Mg-1) (Mg+M) [ (M+l) at  - x] 


1 1/: 
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Table  4 (Cont'd) 


i 


AIRLOAD  EQUATIONS  FOR  CASE  1 (M  < 1,  M > 1) 


Ac  (x,t)  « Ac  (x,t)  + Ac  (x,t)  - Ac  (t) 
P4  p2  p3  pl 


fc-xl1/2 

AcP s**'0  - fir  *(t> 


*(t)  - Ac  (xVI,t) 


X.._  m c + (M  - l)at 

1 1 + M 


Ac  (x)  » 3- 

p6  [1-M2] 


Table  5 


AIRLOAD  EQUATIONS  FOR  CASE  2 (M- < 1,  0 < M < 1) 


Line 


Equation 


(M  + Mg)  at 
(M  + l)at 


x ■ C + (M 

M + 1 


Ac  (X,t) 
pl 


(M+M_)  [ (1+M) at  - x] 


(l+Mg)X 


[(1-Mg)x]l/2  + { (M+Mg) [ (1+M) at-x] } 

M[l-Mg2]  1/2  [(1+M)  (x-  (M+Mg)  at}]  */2 


Ac  (x / 1 ) 
p2 


(M+M  ) [ (1+M) at-x] 


(1+Mg)x 


Mg  ((l-Mg)x]l/z  + { (M+Mg)  [ (1+M)  at  - X]  } 

M £ 1"M  2 ] 1 / 2 [ (1+M) { (M+M„)at  - X } ] / 
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Table  6 


AIRLOAD  EQUATIONS  FOR  CASE  3 (M  < 1,  -1  < M_  < -M) 


Line 


Equation 


c + (M+Mg) at 
c + (M-l)at 


x - c + (M+l)  at 


[M+Mg)  a 


Ac„  (x,t) 
pl 


Z 1 l/2 


cosh 


irM  [l-Mg*  ] 


_1  j (1+Mg)  (C-X) 

I (1-M)  [ (M+Mg) at+c-x] 


Ac„  (x , t ) 
p2 


7TM(1-Mg2]  1/2 


cosh"1  1 - 


( l+M)  (C-X) 


(1-M) ( (M+Mg) at+C-x] 


Ac  (x,t) 
p-> 


TTMd-Mg2]  */2 


cosh 


(M+M  ) at  + c 


(M+Mg) at+c-x 


Table  6 (Cont'd) 


AIRLOAD  EQUATIONS  FOR  CASE  3 (M  < 1,  -1  < M < -M) 


Table  7 

AIRLOAD  EQUATIONS  FOR  CASE  4 (M  < 1,  M < -1) 


Line 


Equation 


- c + (M+Mg) at 
■ c + (M-l)at 


x * c + (M+l) at 

9 


M^+2 (M+M  ) -1 

x * 2 c + (M-l)at 

(M+M  ) (M+l) 


x * c + (M+l) at 


VIII 


(Mg-1) (1-M) 

2 (M+M)  (1+M) 
9 

(Mg-1) (1-M) 
2 (M  +M)  (1+M) 

g 


Ac  (t) 
pi 


4lMgl 

M(M  ! -l! 1/2 


Ac  (X,t) 
p2 


I 2 - 

AcPl(x'fc)  j1  ‘ 7 tan 


(M+M)  (x-at(M-l)-c) 


(Mg+1)  (c-x) 


Table  7 (Cont'd) 

AIRLOAD  EQUATIONS  FOR  CASE  4 (M  < 1,  M < -1) 


Ac  (x,t) 
p3 


Ac  (x,t)  1 - \ tan”1 
p1  ir 


(M+M_) [at (1+M) -x]  + C (M+l) 


(Mg  “ 


2 M 
if  M_  (1+M) 

g 


(M  -1)  13/2 

j"(M+M  ) [at (1+M)  -x]  + C (M+l)  1 

A L 9 * 


Ac  (x,t)  ■ Ac  (x / t } + Ac  (x,t)  - Ac  (x,t) 
P4  P2  P3  Pi 


1C  (x,t)  . 


<J>(t)  * Ac 


VC-  x _ 


XVT  * c + (M-l)at 

1 + M 


Region  6 

Ac  <x,t) 
?6 


where 


¥(t> 


(xVII' fc) 


*VII 


,1/2 


c - 


*VII 


“vil 


(M+l) at  - 


M 

1 - M 


Region  7 


tep^x.U  - acp6  <x'tVXII)  + 


| (yi)  (1-m) | 

^III  * |IhS  2 (Mg+M)  (1+M)  j 


c 

a 


k(x)  ■ Ac  (x)  - Ac  (x,t..TTT) 
Pg  P6  -VIII 


Region  8 


Ac  (x)  3 4 


c-x 

( 1-M2 ) x 


(t"tVIII) 


1/2 


Table  8 


AIRLOAD  EQUATIONS  FOR  CASES  (M  > 1,  M > 1) 


Table  9 


Line 


Equation 


(M+Mg) at 
(M+l) at 
(M-l)at 


Ac  (x,t) 
P1 


Mx-at (M2-l) 


*/* 


(Mz-1) 


(1-M|) 1/2 


cosh 


at (M  M+l)  - M x 


(Mg+M) at-x | 


Ac  (x,t)  * Ac  (x,t) 

P2  Pi 


Ac  (x,t)  * Z 

p3  M ( 1-M2 ) l/2 

9 
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Table  10 

AIRLOAD  EQUATIONS  FOR  CASE  7 (M  > 1,  M < -M) 


Line 


Equation 


I 

II 

III 

IV 


x ■ 0 
X ■ c 

x ■ c + (M+M)at 
9 

M+l 

x m c + (M+l) at 

M+M_ 


x - 51zi_  c + (M-l)at 
M+M_ 


Region  1 

Ac  (t)  - — 
pl 


4M 


M(Mg2-l)  1/2 


Region  2 


Ac  (x,t) 
p2 


4 

1 M ctJ3-l 

'(M+Mg)  (Mx-at(M2-l)  J-(M2-1)c" 

ttM 

fm2-i)  1/2 

(M+Mg)x  J 

M , r(M+M  ) { (1+HM  ) at-M  x/+ (1+MM  ) C ] | 

— 2 2 — m cos"  2 2 2 2 , 

(M  -1)  ' L (M+Mg)  {c-X+ (M+Mg)  at } J; 


Region  3 
Ac  (t) 


P3  (M2  -1) */2 


51 


(24) 
Vt 

+ — , leading  edge  intercepted  first 

(25) 
Vt 

+ — , trailing  edge  intercepted  first 

12 

The  dynamic  pressure  q “ and  ambient  pressure  pa 

used  to  calculate  the  local  pressure  on  the  airfoil  are  based  on  the 

preblast  values  of  p^,  pa,  and  until  Acp  becomes  nonzero.  After 

that  time,  the  blast  values  of  p , p and  V at  the  leading  edge  of 

00  00  £ 

the  airfoil  are  used. 


(1)  Subsonic  Flight  Speed  (M  < 1) 


(a)  Case  1:  M > 1 

S 


The  equations  for  Acp  in  Case  1 were  derived  by  Smiley 
and  Krasnoff  (ref.  19).  The  shock  wave  first  intercepts  the  wing  at 
leading  edge  and  sweeps  over  the  airfoil  at  a speed  greater  than  sound 
as  shown  in  figure  12.  Region  1 lies  between  the  shock  wave  III  and 
the  rearward  emanating  wave  IV.  These  waves,  in  turn,  reflect  waves  V 
and  VI  when  they  reach  the  trailing  edge,  forming  Regions  2,  3,  and  4. 


Extensive  calculations  are  required  to  determine  Ac^ 

in  Region  5,  yet  the  distribution  of  Acq  is  rather  simple  in  form, 

blending  from  the  value  at  reflected  wave  VI  to  the  zero  value  at  the 

trailing  edge.  The  loading  in  this  region  is  relatively  unimportant, 

so  the  steady-state  distribution  for  Ac^  is  obtained  by  sea" xng  to 

match  Ac  at  reflected  wave  VI.  When  reflected  wave  VI  reaches  the 
P 

leading  edge,  a steady-state  pressure  distribution  exists  over  the 
entire  airfoil. 
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(b»  Case  2:  0 < M <1 

g 

The  solution  for  Acp  for  Case  2 was  derived  by  Ruetenlk 
(ref.  14)  for  Regions  1 and  2.  No  solution  has  been  obtained  for 
Region  3 due  to  the  complexity  of  the  interaction  between  the  singular 
distributions  of  loading  in  Regions  1 and  2 with  the  wake.  The  distri- 
bution of  40^  for  Region  3 is  obtained  by  using  the  steady-state  dis- 
tribution for  Ac  and  matching  with  Ac  at  reflected  wave  V.  The  match- 
P P 

ing  is  accomplished  by  a linear  stretching  of  the  coordinates,  based 
on  distance  from  the  trailing  edge,  to  match  the  loading  with  the 
loading  along  wave  V.  The  specific  equations  used  are 


x 


c 


. Jit). 

C - xv(t) 


(c  - x) 


with  x(t)  defined  by 


(26) 


(27) 


(28) 


where  Ac 

Pl. 

as  appropriate,  and  x^Ct)  is  the  chordwlse  position  of  wave  V.  The 
steady-state  distribution  of  Ac^  is  reached  when  wave  V intercepts  the 
leading  edge  of  the  airfoil.  This  concurs  with  the  time  required  to 
reach  steady-state  as  indicated  by  equation  (59),  reference  21.  The 
steady-state  distribution  for  Acp  applies  after  wave  V intercepts  the 
leading  edge. 


(Xy,t)  is  the  value  of  Ac  along  wave  V for  Region  1 or  2, 


(c)  Case  3:  -1  < M ^ -M 

g 

In  'Jase  3,  the  shock  arrives  from  cha  trailing  edge, 
traveling  forward  faster  than  the  airfoil  Mach  nuaber  but  slower  than 
the  speed  of  sound.  The  wave  diagraa  is  shown  in  figure  14. 

The  equations  for  for  Regions  1,  2 and  3 were  derived 
using  the  transformation  from  the  two-dimensional  supersonic  steady- 
state  wing  to  the  transient  one-dimensional  airfoil  along  with  the 
equations  for  the  loading  given  in  Section  A13,  reference  11.  The  load- 
ing in  Region  4 is  nearly  the  same  as  in  Region  2,  except  for  a slight 
perturbation  from  Region  3.  Because  the  loading  in  Region  2 is  small, 
the  ACp  value  for  Region  4 is  obtained  using  the  equation  for  Region  2. 

The  steady-state  loading  distribution  is  used  for  Region  5, 
scaling  the  distribution  asymptotically,  following  equation  (18), 
reference  21.  All  values  of  Acp  in  Region  5 are  scaled  by  the  factor 
*(t). 


* - l - 2-T7  (29) 

T - i 

c M + M 

8 

The  time  factor  HO  varies  from  0.5  when 
leading  edge,  to  a maximum  of  1.0. 


(30) 

the  blast  shock  reaches  the 


(d)  Case  4: 


< 


-1 


In  this  situation,  the  blast  shock  overtakes  the  airfoil 
from  the  rear  at  a speed  greater  than  the  speed  of  sound.  The  equation 
for  ACp  for  Regions  1 through  4 are  taken  from  Smiley  and  Krasnorf 
(ref.  19).  In  Regions  5 and  6,  the  equation  for  the  steady-state  dis- 
tribution of  Ac  is  matched  to  Region  4. 


(2)  Supersonic  Flight  Speed  (M  > 1) 


The  equations  for  all  regions  in  Cases  5,  6,  and  7 were 
obtained  from  table  3 of  reference  20. 

3.2  LOADING  ON  FUSELAGE  SURFACES:  PREBLAST  AND  BLAST  ENCOUNTER 

Since  the  fuselage  shape  is  not  known  a priori,  it  is  difficult 
to  set  up  a practical,  yet  rational,  procedure  for  computing  the  time 
variation  of  the  pressure  at  an  arbitrary  point  on  the  actual  fuselage. 
The  devised  procedure  is  an  attempt  to  estimate  the  pressures,  admit- 
ting the  fact  that  much  of  it  is  based  on  intuition  and  Involves  many 
simplifications  (which  have  been  adopted  in  related  situations,  e.g., 
see  Norris  and  Hansen,  reference  22)  which  may  not  be  fully  justifiable 
for  certain  configurations. 

If  wind  tunnel  data  were  available  for  the  particular  fuselage, 
such  data  could  be  used.  These  data  would  have  to  include  pressure 
coefficients  at  a sufficient  number  of  points  on  the  fuselage  for  a 
range  of  Mach  numbers  and  Reynolds  numbers.  The/  would  also  have  to 
cover  a wide  range  of  angles  of  attack  and  angles  of  yaw,  and  combina- 
tions thereof.  Unfortunately,  such  extensive  data  are  never  obtained 
for  actual  configurations.  In  addition,  the  variety  of  configurations 
encountered  prevents  the  use  of  one  set  of  data  for  all  configurations. 

Other  theoretical  means  are  available.  Methods  have  been  devised 
to  determine  pressures  on  bodies  at  suosonic  and  supersonic  speeds 
for  certain  types  of  bodies.  Methods  for  subsonic  and  supersonic 
cylindrical  bodies  (e.g.,  refs.  23-26)  are  not  only  limited  but  also 
extremely  laborious.  They  are  rejected  here  on  practical  grounds 
recognizing,  in  addition,  that  other  aspects  of  the  problem  can  only 
be  treated  in  very  approximate  fashion.  The  Munk-Jones  slender-body 
theory  or  the  linearized  supersonic  theory  for  pointed  bodies  of  revo- 
lution will  be  of  some  help.  The  Munk-Jones  theory  gives  running  loads 


per  unit  body  length,  independent  of  Mach  number,  end  considers  only 
chose  additional  pressures  due  to  angle  of  attack.  It  has  certain 
complications  if  the  cross  section  is  not  circular.  Both  of  the  slender- 
body  theories  (Munk-Jones  and  llnsarizad  supersonic  theory  for  pointed 
bodies  of  revolution)  are  Mach  number  independent  for  the  pressures 
induced  by  angle  of  attack.  The  pressures  at  zero  angles  of  attack  can 
be  obtained  for  the  supersonic  case;  however,  they  are  Mach  number 
dependent,  and  cannot  be  easily  introduced. 

The  point  is  that  neither  sufficient  data  nor  well-established 
practical  methods  exist  to  perform  even  the  steady-state  calculations 
of  pressures  called  for  in  the  present  problem  for  the  preblast  and 
post-diffractive  phases  of  the  pressure-time  histories.  The  rough 
approximations  which  are  described  below  appear  to  offer  the  best 
solution. 

Consider  the  two  bodies  in  figure  19,  where  the  actual  fuselage, 

(a),  and  an  equivalent  body  of  revolution,  (b),  have  the  same  cross- 
sectional  area  distributions  along  their  lengths.  For  purposes  of 
equating  the  pressures  on  the  two  bodies,  points  on  the  two  bodies  may 
be  related  according  to  either  of  the  two  following  schemes 


(a) 


yb  yb 

9 — 9,  i.e.,  ~ ~ T- 

7b  S 


XB  ' '*B 


This  implies  chat  cha  points  which  have  tha  same  slope  of  the  contour 
in  the  plane  section  B-B  are  related. 

Scheme  (b)  was  chosen,  since  it  results  in  the  proper  inclination 
to  the  relative  flow.  Selection  of  the  point  (x^,  y^ , z^)  where  the 
pressure  versus  time  is  desired  permits  unique  determination  of  the 
related  point  (x^,  , zfc)  on  the  equivalent  body  where  the  pressures 

will  be  obtained. 

The  replacement  of  the  actual  body  by  a body  of  revolution  repre- 
sents a very  substantial  simplification.  If  the  direction  of  the  blast 
is  arbitrary,  the  blast  will  produce  an  incremental  angle  of  yaw  (a^) 
in  addition  to  an  Incremental  angle  of  attack.  However,  these  two  can 
be  combined  with  the  initial  angle  of  attack  (oq)  in  the  preblast  pitch 
plane  (z^  - x^)  to  form  a single  angle  of  attack  in  a plane  different 
from  the  pitch  or  yaw  plane. 

Consider  the  equivalent  body  and  the  cross  section  for  a point 
(x2b»  y2b • z2b^  where  the  pressure  is  desired,  as  shown  in  figure  20. 
Note  that  the  body-coordinate  system  is  left-handed  with  the  origin 
at  any  point  on  the  x^-axis  of  the  body. 

The  objective  is  to  find  the  pressure  at  point  (x^,  y2b’  Z2b^ 
before  and  after  the  blast  encounter.  The  pressure-time  history  will 
be  assumed  to  vary  in  the  fashion  shown  in  figure  21.  This  is  similar 
to  the  pressure-time  history  presented  by  Norris  and  Hansen  (ref.  22) 
for  stationary  cylinders  encountering  blasts. 

The  pressure  variation  spans  the  following  three  time  ranges  of 
interest : 


(a)  Period  prior  to  blast  arrival,  t <_  t where  t denotes 
the  time  of  shock  arrival  at  (x,,^.  ^2b’  Z2b^ 

(b)  The  diffraction  period,  t < t < t.,  where  t.«t  +t 

a d a c a 

denotes  the  time  when  the  flow  returns  to  essentially 
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"steady-state"  condition*  following  blast  «ncount*r  and 


1 j • 


t is  tha  daaring  time.  An  abrupt  increasa  in  pressura 
may  occur  at  time  t+  (or  t + e,  e ■*  0). 

A A 


(c)  Tha  post-diffraction  phasa  or  tha  drag  phasa,  t >_  t^. 
Ranga  (a)  is  truly  a staady-stats  period.  Range  (c)  can  be 


treated  as  steady-stata  basad  on  instantanaous  fraa-straaa  conditions 
dictated  by  blast  conditions  at  the  point  in  question.  To  construct 
pressure  versus  time  plots  like  that  of  figure  21,  the  following  prob- 
lems need  to  be  resolved: 


1. 


Determination  of  tha  steady-state  pressure  at  the  point 
of  interest  based  on  preblast  conditions. 


2. 

3. 

4. 

5. 


The  shock  arrival  time,  t#,  for  the  point  of  Interest. 


Tha  length  of  tha  diffractive  period  t^. 


Tha  pressure  at  time  t+,  p(t+). 

A tt 


Tha  pressure  at  time  t^,  p(t^),  based  on  Instantaneous 


"gust"  and  "free-stream"  conditions  at  td» 


6.  Pressures  during  the  diffractive  period,  assuming  linear 


variation  with  time  between  pressures  p(t  ) and  p(t,). 

A G 


7. 


Pressures  at  selected  times  t > t^,  p(t),  based  on 


instantaneous  "gust"  and  "free-stream"  conditions  at  t. 


Each  of  these  will  now  be  discussed  individually.  First,  however, 
some  preparatory  work  is  needed.  Divide  the  body  into  three  parts, 
as  shown  in  figure  22.  The  definitions  are  arbitrarily  made  so  that 


is  the  station  at  which  (dR/dx^) 


-0.15  and  x is  the  station  at 
S 

which  (dR/dx^)  ■ 0.15,  where  R is  the  radius.  Pressures  on  Part  (c) 
are  not  expected  to  be  large,  especially  for  large  forward  speeds, 


V ; therefore.  Part  (c)  can  be  treated  in  a very  approximate  fashion  as 
o 


a portion  of  Part  (b)  with  (dR/dx^)  > 0. 
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Definition  of  Parts  of  Equivalent  Body 


3 

1 


The  steady-state  pressure  coefficient  (C^),  due  to  angle  of  attack 
for  the  quasi-cylindrical  section,  is  a function  of  position  (x^  and 
6)  and  angle  of  attack  (a)*  according  to  slender-body  theory. 


/ 2 dR  . , ...  . 2 dR 

4acos  a 7 — cosQ  + a (A)  sin  a S.0 


C (V9,a) 


(31) 


,a(A)  sin  a 


dR 

dxb 


> 0 


.7-  9 a > 0 


9 


a < 0 


(32) 


l-4sin  A 0 < A < tt/2 

a(A)-<  -10.278  + 4.6333A  tt/2  < A < 2tv/3 

-0.574  2tt/3  < A < it 


(33) 


The  term  4a  cos  a(dR/dx^)  cos6  comes  from  slender-body  theory  (see 

Liepmann  and  Roshko,  reference  25,  page  245,  with  0}  omitted  and 
2 

a cos  a added  to  reduce  the  total  dynamic  pressure  to  dynamic  pressure 

2 

based  on  axial  velocity.  The  added  cos  a factor  also  reduces  the 
(dR/dx^)  contribution  for  large  angles  of  attack.  The  a (A)  sinfca 


* 

The  angle  of  attack  is  defined  as  the  angle  between  the  negative  of 
the  relative  wind  and  the  body  axis  x,  . 
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term  gives  Che  cross-flow  drag  in  accordance  with  Allen  and  Perkins 
(reference  26).  The  constants  in  the  formulas  for  the  circumferential 
pressure  variation,  a(A),  are  chosen  so  that  the  circumferential  pres- 
sure variations  result  in  a value  of  0.35  for  the  sectional  coefficient, 
C^,  a reasonable  value  for  supercritical  Reynolds  numbers.  The  circum- 
ferential pressure  variation,  shown  in  figure  23,  matches  a typical 
pressure  variation  along  the  periphery  of  a cylinder  engulfed  in  a 
supercritical  Reynolds  number  flow. 

Equations  (31)  describe  the  steady-state  pressure  coefficient  for 
the  quaai-cylindrical  and  tail  sections  of  the  body,  but  are  not  appro- 
priate for  the  blunt  nose  section.  The  nose  section  pressures  will  be 
approximated  by  using  the  pressure  coefficient  at  x^x^,  which  will  be 
taken  as  the  stagnation  point.  This  pressure  coefficient  is  a function 
of  Mach  number,  M,  only,  and  will  be  denoted  by  C . 


(l+M2/5)7/2-l 


M < 1 


0 . 7M 


C (M)  - 
Po 


(34) 


7/2 


5/2 


(*)'  fe)  - - 


0.7M 


Consider  now  the  seven  problems  listed  earlier.  The  first  is  that 
of  determining  the  preblast  pressure.  For  the  quasi-cylindrical  and 


tail  sections,  C (t  < t ) is  given  by  equation  (31),  using  the  pre- 
P ~ a 


blast  angle  of  attack. 


The  determination  of  the  pressure  coefficient  for  the  nose  section 
is  more  complicated.  At  the  nose  point,  which  is  the  stagnation  point 


for  a ■ 0,  C is  applicable  exactly  at  a ■ 0,  and  approximately  if  a 
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is  finite.  The  pressure  coefficient  for  a quas '.-cylindrical  section, 
Cp(xd,9,°0’  is  defined  by  equation  (31)  where  9 is  the  same  as  that 
for  point  (xb,  yfe,  zb). 


(35) 


(36) 


The  transition  from  the  stagnation  point  pressure  coefficient  to  the 
pressure  coefficient  at  x^  for  t < t can  be  represented  by 


I sin6  - sin6 

V W6^  ’ [Cp0  ' Vve>c,)1  [-  : 


+ C (x  ,e,a) 
P d 


(37) 


The  form  of  equation  (37)  is  patterned  after  the  pressure  distribution 

2 

given  by  Newtonian  flow,  that  is,  proportional  to  sin  5.  The  variation 

in  pressure  coefficient  given  by  equation  (37)  roughly  approximates  this 

distribution  because  C is  substantially  larger  than  C .and  sin5,  is 

PQ  P d 

small  in  comparison  with  unity. 

The  pressure  coefficients  defined  above  are  used  to  obtain  the 
pressure  in  accordance  with  the  equation 


P 


p + q c 

• n r 


(38) 
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where  p^  and  q^  £ Che  free-stream  pressure  and  dynamic  pressure. 

Equation  (38)  is  appropriate  for  transforming  the  pressure  coefficient 
into  pressure  tor  both  preblast  and  blast  conditions,  provided 
that  p and  q describe  the  local  ambient  conditions  in  front  of  or  at 
some  point  behind  the  shock,  as  appropriate.  One  note  of  caution  is 
required  in  applying  equation  (38):  The  pressure  calculated  by  equation 
(38)  may  be  negative,  which  is,  of  course,  a physical  impossibility. 

If  this  happens,  the  pressure  should  be  set  equal  to  zero. 

Next,  the  determination  of  the  time  of  shock  arrival,  t , and 

the  clearing  time,  t , must  be  considered.  The  blast  characteristics 

c 

routine  used  in  the  NOVA  program  returns  the  shock  position  for  a given 
time.  It  is  a simple  matter,  then  to  proceed  in  small  timewise  steps, 
calculating  the  position  of  the  desired  point  on  the  body  relative  to 
the  burst  point  at  each  step,  and  comparing  this  with  the  shock  position 
at  the  time  associated  with  that  step.  In  this  way,  the  time  of  arri- 
val, t , can  be  identified  as  closely  as  desired, 
a 

Similarly,  the  time  at  which  the  sho'-k  first  touches  the  fuselage 
can  be  identified.  The  time  at  which  the  fuselage  is  intercepted  t^,  ia 
needed  to  define  the  clearing  time,  t^.  Determination  of  t^  is  some- 
what more  cumbersome  than  the  determination  of  t , because  a number 

a 

of  points  on  the  fuselage  must  be  considered  to  identify  the  time  at 

which  the  shock  first  intercepts  the  fuselage.  3asically,  however,  t^ 

is  determined  in  the  same  way  as  t . 

a 

The  clearing  time  is  an  artificial  parameter  wnicn  has  been  devised 
to  present  the  actual  pressure  versus  time  in  the  approximate  fashion 
of  figure  21.  It  may  be  thought  of  as  a rough  measure  of  the  time 
required  for  the  flow  to  return  to  essentially  quasi-steady  conditio'-'' 
following  shock  arrival.  Attempts  to  define  t experimentally  are 
rather  limited  and  are  confined  primarily  to  tests  on  stationary  cylinders 
reference  27). 


There  is  sone  degree  of  arbitrariness  in  defining  t because  (a}  the 

c 

diffractive  pressure  variations  are  not  generally  linear,  and  ib/  the 
start  of  the  quasi-steady  period  aav  be  very  difficult  tc  deteraine. 

Both  of  these  points  are  strongiy  influenced  by  the  location  of  the 
point,  the  shape  of  the  body,  and  the  shock  orientation  and  strength. 

In  what  follows,  then,  an  atteapt  is  aade  to  give  a definition  of  : 

c 

which  will  result  in  reasonable  representations  of  the  diffractive 
phase  pressures.  There  is  obviously  no  single  approach  which  ccu-a 
handle  accurately  all  poss)w’_  :ases  that  aav  arise. 

To  begin  with,  co>'.  i-er  t.ie  ase  of  a wedge  flying  supersonically 

prior  to  a head-on  bla'  en„uur.tar.  During  the  blast  traversal,  the 

snock  pattern  has  been  noted  to  ne  such  tnat  a region  behind  t.ne  shock 

of  approximately  a^U-t^)  in  width  )a^  “ speed  of  sound  behind  the 

blast  shock)  is  still  in  a diffractive  pnase  (re:.  13).  Figure  2- 

illustrates  this  point.  Thus,  one  =av  consider  for  the  above  situation 

that  a "diffraction  wave"  starts  traveling  back  from  the  shock  wave 

at  tine  t..  This  wave  travels  at  the  speed  of  souna,  a.^^,  relative  to 

the  shock.  At  some  tine  t > t.,  then,  the  diffraction  wave  is  at  a 

distance  a.  (t-C.)  behind  the  shock.  Since  the  position  of  the  diffrac- 
ts i 

tion  wave  can  be  referenced  to  the  position  of  the  shock  front,  the 
procedure  described  above  can  be  used  to  deteraine  the  tine  at  whicn 
the  diffraction  wave  reaches  the  desired  point,  t . The  clearing  tine, 

t , is  then  siarplv  the  difference  between  t.  and  t . 

c da 

If  this  analytical  cerimtion  for  t is  adopted  for  bodies  flying 

c 

at  different  speeds,  and  encountering  blasts  f reus  arbitrary  orienta- 
tions, some  problems  arise.  First  consider  the  case  of  a cylinder 
which  is  intercepted  side-on  by  a shock  wave.  The  above  definition  of 

t results  in  a clearing  tine  of  zero  for  the  extreme  windward  points, 
c 

since  t.  and  t are  identical  for  such  points,  and  therefore,  the 
1 a 

diffraction  wave  also  arrives  at  t (t.  ■ t ).  Tests  on  scaled  missile 

ad  a 

models  in  the  DASACON  shock  tube  (ref.  27)  however,  indicate  that  a 
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clearing  time  of  about  (ZR/a^)  is  appropriate,  particulai 'y  for  points 
on  the  windward  side  which  experience  the  highest  pressures.  Accord- 
ingly, the  adopted  procedure  is  to  determine  t^  as  described  earlier, 
and  to  use  the  larger  of  that  value  or  (2R/a^s>. 

The  second  problem  which  arises  may  be  illustrated  by  the  case 
when  the  blast  arrives  directly  from  the  back  and  the  vehicle  velocity 
nears  the  value  Vg  - a^g,  where  Vg  is  the  velocity  of  the  shock  front. 

or  such  cases,  the  clearing  time  ( t > , as  defined  above,  becomes 
unrealistically  high.  No  experimental  results  are  available  to  establish 
an  upper  limit  on  t in  a manner  such  as  that  used  to  establish  the 
lower  limit.  The  value  of  (6R/a^g)  has  arbitrarily  been  chosen  as  the 
upper  limit  of  t in  the  program. 

To  determine  the  pressure  at  time  t*,  it  is  necessary  to  define 
the  orientatioa  of  the  shock  front  relative  to  the  fuselage  local  sur- 
face normal.  When  t is  determined,  as  described  previously,  the  shock 

a 

orientation  in  apace  is  automatically  available.  The  unit  vector 
normal  to  the  shock,  in  the  direction  of  shock  propagation,  will  be 
designated  by  ng.  The  unit  vector  normal  to  the  body  and  directed 
inward  is  given  by 


"b 


Cl  — 

D axb 


■ Vlnc 


k.  cose) 
o 


(39) 


where  i^  jfe, 
tions,  respec 
is  related  to 
axis,  so  that 


and  k^  are  the  unit  vectors  in  the  x^ , , and  z^ 

tivelv.  Defining  a new  coordinate  system,  x,  v,  z 

the  x,  , v.  , z.  svstem  bv  a rotation  —a  about  the 
b b d • • o 

the  x axis  is  horizontal. 


airec- 

whicn 

:/b 


■ I 1 (^—  co»a^  + co*6sina^)  - j sin  ? 


r ,dR 

k (- — siaa 
dxb 


cosdcosa  ) 
o 


VlH%) 


where  i,  j,  end  k are  the  unit  vectors  in  the  x,  y,  and  z directions, 
respectively. 

If  the  components  of  n (the  unit  vector  normal  to  the  shock  in 

s 

the  direction  of  shock-front  propagation)  in  the  x,  y,  z directions  are 
n , n , n , respectively,  then  the  shoe*  plane  and  tne  plane  tangent 

* y * 

to  rhe  fuselage  body  sake  an  angle 


- cos'  . ns) 


• :r.Ta=ss  ■ n — cosa  ♦ cosrsina  ) 

a - x d\  0 


-n  sine  + n (— — sina  -cosccosa  ) V 0 * i , < 
y z dx^  o o / - sb  - 


If  the  vehicle  has  no  component  of  its  velocity  normal  to  the  surface, 
then  the  theory  of  plane  shock  reflection  from  a stationary  plane  sur- 
face can  be  used.  This  would  be  nearly  true  for  regions  where  dR/dx.^ 
is  small,  i.e.,  in  quasi-cylindrical  portions.  Bor  high  dR/dx^ 
regions,  i.e.,  the  nose  section,  this  would  not  be  strictly  applicable. 

In  fact,  with  "blunt"  bodies  in  supersonic  flow,  there  is  a shock-shock 
interaction  problem  which  is  not  easily  Included. 
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la  ch*  absence  of  a better  practical  procedure  for  such  regions, 

the  reflection  theory  of  a plane  shock  on  a plane  stationary  wall  will 

be  used.  Norris  and  Hansen  (ref.  22)  give  Ap/Ap^  as  a function  of 

9 , where  Ap  is  the  increaental  pressure  on  the  body  and  Apg  is  the 

shock  overpressure.  Ap/Apg  versus  9^  is  dependent  on  the  pressure 

ratio  across  the  shock.  The  curves  of  Ap/Ap^  versus  9 ^ exhibit  buaps 

in  the  neighborhood  of  transition  froa  a regular  reflection  to  a Mach 

reflection.  These  buaps  have  not  been  observed  experiaentally  on  'ircular 

cylinders.  father  than  use  the  curves  in  Norris  and  Hansen  directly,  the 

approximation  depicted  in  figure  25  is  used.  The  curve  shows  constant 

Ap/-p  for  0 ■ 7 . - and  linear  variations  of  ip/ip  with  •?  , between 

* “ so  ^ s sb 

and  ’12,  and  between  • / 2 and  *. 

In  figure  25,  the  parameter  r^  is  the  reflection  factor  for  normal 
incidence,  ard  is  given  by 


P,  * 4/7  ^ 

P.  * 1/'  Apc 


Thus,  referring  to  figure  23, 


0 < 9 < t / 4 

— sb  — 


r + u ■ r ) ( 

p p \ 7/4 


1.5  - 


( y74  | 1 < , < t/2 

\ 7/4  / 4 sb  — 


— 9 <7 

2 sb  - 


The  pressure  at  t is  then  determined  from 


iNCOEMCE  XN6LE, 


Figure  25.  Variation  of  Reflection  Factor  w 
Incidence  Angle 


P(t*>  - >>u;>  ♦ ip. 


where  p(t”)  is  Che  pressure  just  before  shock  arrival,  or  che  original 
ste.idy-scate  pressure  froei  equations  (31)  through  (38). 

Next,  consider  che  pressure  ac  che  end  of  che  diffraction  period, 

p(t  ).  Proa  che  blast  routines  and  che  range  and  oriencatlon  of  che 

point  in  question  relative  to  che  burse  poinc,  che  aabient  density, 

-^(t^),  pressure,  P.U^),  and  aaterial  velocity  coaponents,  v^t^), 

w (t.),  and  w (t,).  are  known.  The  x,  y,  z coordinate  systea  is  that 
yd  z d 

described  in  connection  with  equation  (40).  The  total  velocity  is, 
chen. 


vv  *vv<.  - “«!vr  * lyv  * 'w 


To  find  che  angle  of  attack  of  che  section  containing  che  point  in 
question,  che  total  velocity  of  the  air  relative  to  the  body  is  resolved 
into  coaponents  in  che  x^,  y , z^,  body  coordinate  system. 


V„  - - (V  - w ) cosa  + w stna 

T ox  o z o 


VT  “ Wy 
yb 


V„  * (V  - w ) sina  + w cosa 
T ox  o z o 

2b 


The  flow  is  then  inclined  at  an  angle  y from  the  z^  axis,  where 
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An  effective  9 nay  be  defined  ee  the  angle  between  the  flow  and  the 
point  in  queation.  Calling  thia  angle  9, 

9 i 9 - y (50) 

To  maintain  9 in  the  interval  -f  - f,  if  9 > t,  subtract  2*.  and  if 
9 ' add  2*.  Finally,  the  angle  of  attack  is  given  by 


i • tan 


-1 


V,  ain>  ♦ V cost 

i I 


T 


a 


< t/2 


(51) 


Equations  (31)  through  (38)  are  used  with  9 replacing  9 to  find 
p(t.).  Note  that,  if  a nose  point  is  involved,  there  is  an  implicit 
assumption  that  * 0.  Also,  for  a nose  point,  che  flow  parameters 

*b 

to  be  employpH  for  the  stagnation  point,  including  the  Mach  number  in 
equation  (31),  and  for  . are  those  corresponding  to  the  point 

in  question. 

Pressures  during  the  diffraction  period  are  found  by  linear  inter- 
polation. 


p(t) 


P(0  + 


t-t 


t 

c 


t 

a 


- Cd 


(52) 
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For  daes  greater  chan  t^,  the  pressure  at  a point  is  found  in 

precisely  Che  saae  way  as  p(t^);  chat  is,  equations  (45)  through  (51) 

are  used,  with  t > t,  replacing  t.,  and  the  parameters  so  defined  are 

d d 

used  with  equations  (31)  through  (38)  to  obtain  the  pressure. 


The  foregoing  procedures  serve  to  define  coxpletely  the  pressure- 
ciae  history  for  a point  on  a fuselage  with  a straight  centerline,  due 
to  a single  shock.  Modifications  necessary  to  accomodate  a nonscraight 
centerline  and  a double  shock  will  now  be  described. 


The  computer  program  permits  the  centerline  of  the  fuselage  to  be 

described  by  two  straight  lines,  as  shown  in  figure  26.  The  bend 

occurs  at  ■ x^ , which  is  assumed  to  be  aft  of  the  transition  point 

from  the  nose  section  to  the  quasl-cvlindrical  section,  x..  All  of 

d 

Che  preceding  procedures  remain  valid,  as  long  as  the  initial  angle  of 
attack  is  defined  bv 


% 


BF 


a - 
o 


(53) 


*b 


*BF 


If  the  aircraft  is  intercepted  by  the  blast  wave  at  a point  just 

above  the  triple  point,  a second  shock  will  follow  closely  behind  the 

first.  The  pressure-time  history  for  this  case  is  shown  in  figure  27. 

Subscripts  1 and  2 are  used  to  identify  the  first  and  second  shocks, 

respectively.  Prior  to  the  time  t , the  pressure  is,  of  course,  given 

a2 

by  the  procedures  already  described.  The  jump  in  pressure  at  time  t 

a2 

is  (Ap/Ap  ) Ap  . The  shock  overpressure  Ap  , for  the  second  shock 
S2  a2  S2 
is  obtained  from  the  blast  routine  and  (Ap/Ap  ) is  obtained  from 

S2 
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Figure  26.  Geometry  of  Fuselage  with  Bent  Centerline 


PRESSURE 


equations  (42)  and  (43)  using  Che  ambient  pressure  corresponding  to  time 

t for  p . The  pressures  at  t and  later  are  given  by  procedures 
a2  " a2 

already  described,  using  flow  parameters  from  the  blast  routine  for  the 

doubly  shocked  air.  The  only  remaining  point  requiring  definition  is 

p(t  ) for  t < t . 
dl  a2  dl 


Cd  * Ca 

p(t  ) . p <C  ) + (p(t*  > - P<t"  ))  - - ■ 1 . - ? 

1 1 2 2 d,  a. 


(p(t  ) 
a2 


■P(ta..))  * (p(td  > - Pl^d  >> 


c < t , 
a2  dl 


(54) 


The  quant -ty  p.  (t  ) in  equation  (54)  is  the  pressure  that  would  exist 

1 

at  time  t in  the  absence  of  the  second  shock  and  p,(t.  ) is  to  be 
\ 1 d2 
interpreted  similarly.  Linear  variations  in  pressure  are  assumed  from 

t to  t.  and  from  t.  to  t.  . 
a2  d1  dx  d2 


SECTION  IV 
STRUCTURAL  RESPONSE  TECHNIQUES 


The  structural  elements  of  an  aircraft  vehicle  susceptible  to 
overpressure  damage  consist  of  the  following: 

1.  Stringers  on  the  fuselage  or  on  lifting  surfaces 

2.  Longerons  or  frames  on  the  fuselage 

3.  Single-layered  or  honeycomb  panels  on  the  fuselage  or  on 
lifting  surfaces 

4.  Acrylic,  glass,  or  plexiglass  canopies 

5.  Multi layered  radomes 

6.  Rib  webs  on  lifting  surfaces 

Two  structural  response  codes  were  developed  for  predicting  the 
elastic  and  inelastic  response  of  these  structural  elements  to  the 
overpressure  loading  associated  with  the  blast  wave.  The  first  of  these 
response  codes  is  called  QEPROB  for  Dynamic  Elastic  Plastic  Response  of 
Beams  and  was  developed  for  use  with  stringers,  longerons,  frames,  ribs, 
and  also  for  conical  or  cylindrical  shaped  radomes  which  may  be  repre- 
sented bv  a ring.  The  second  of  these  response  codes  is  called  DEPROP 
for  Dynamic  Elastic  Plastic  Response  of  Panels  and  was  developed  for  use 
with  single-layered  and  honeycomb  panels,  canopies,  and  certain  radomes. 

Provision  is  made  in  each  code  to  determine  the  static  preblast 
solution  due  to  internal  pressurization  and  eerodynamic  loads  on  the 
structural  element.  This  steady-state  solution  then  provides  the 
initial  conditions  for  the  dynamic  response  associated  with  a time- 
dependent  blast  wave.  The  formulation  for  each  of  the  two  response 


codes  is  described  in  detail  in  the  remainder  of  this  section. 


4.1  DYNAMIC  ELASTIC-PLASTIC  RESPONSE  OF  BEAMS  - DEPROB 


The  structural  response  program  DEPROB  calculates  the  static  and 
dynamic,  elastic  or  elastic-plastic,  response  of  aircraft  elements  which 
can  be  modeled  as  beam  or  ring  elements.  These  elements  can  have 
arbitrary  spanwlse  shape,  can  have  arbitrary  cross  section  involving 
different  materials  and  varying  spanwlse  along  the  beam,  can  have  any 
combination  of  boundary  conditions — clamped,  simply  supported,  or  free — 
and  can  respond  to  a transient  pressure  function  which  varies  with  both 
time  and  spanwlse  position. 

The  basis  for  the  code  is  a finite-difference  method  developed  by 
MIT  (refs.  29  and  30)  to  predict  the  deformations  of  a circular  ring  to 
impulsive  and  transient  loadings.  DEPROB  represents  considerable  modi- 
fication and  extension  of  the  capabilities  of  that  original  effort.  The 
remainder  of  this  section  will  describe  in  detail  the  DEPROB  code.  As 
partial  verification,  the  final  portion  presents  comparisons  of  the 
dynamic  response  of  two  clamped  beams  tested  experimentally  and  analyzed 
by  DEPROB  and  by  the  MIT  code. 

4.1.1  Basic  Theory 

The  finite  difference  technique  applied  to  two-dimensional 
structures  in  DEPROB  assumes  a spanwlse  model  consisting  of  a series  of 
discrete  masses  interconnected  by  straight,  weightless  bars,  as  indi- 
cated in  figure  28.  Beams  with  variable  geometrical  cross  section  in 
the  spanwlse  direction  are  reduced  to  a series  of  links,  each  completely 
uniform.  Each  bar  and  mass  then  has  its  own  material  properties.  Beam 
response  to  externally  applied  forces  is  computed  at  each  mass  point. 

The  actual  cross  section  is  idealized  by  introducing  a set  of 
discrete  points  called  flanges.  These  flanges  have  normal  stress 
distributions,  and  are  interconnected  by  material  of  infinite  shear 
rigidity  as  indicated  in  figure  29.  Flange  elongations  are  governed  by 
the  Bernoulli  assumption  that  plane  sections  remain  plane. 
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The  stress-strain  curve  of  each  layer  of  material  is  modeled 
by  a series  of  straight  lines,  as  shown  in  figure  30.  Enough  straight 
line  segments  will,  in  the  limit,  describe  any  curve.  The  curves  may 
differ  in  tension  or  compression,  except  for  the  elastic  slopes  (segment  1), 
which  must  be  the  same.  However,  the  yield  stress,  defined  as  the  break 
point  stress  at  the  end  of  the  first  segment,  may  differ  in  tension  and 
compression. 

A more  complete  description  of  the  theoretical  development  of 
DEPROB  is  now  presented  in  subsections  4.1.2  to  4.1.11. 

4.1.2  Equations  of  Motion 

The  governing  equations  of  dynamic  equilibrium  take  the  form 


(N  cos  9)  - ~ (Q  sin  9)  + Fy  - mV 


- 0 


(N  sin  9)  + ~ (Q  cos  8)  + F - 
9s  9s  z 


and,  ignoring  rotary  inertia, 


Referring  to  figure  28,  the  corresponding  finite  difference  equations 
are,  for  the  ith  mass  point, 


Ni+icos0i+rNicoserQi+isin0i+i + Qisin0i 


Asi  + Asi+ll 

+ Fyi  l — 2 r Vi " 0 


(55) 
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Nl+lsin9i+l“Ni3in9i  + Qi+lcos9i+l 

/Asi  + As^  A 

-QiCO90i + FZi  \ — r-^r  miwi  - 0 

M,  - M.  . - Q,  As.  - 0 
i i-1  1 


where,  from  geometrical  considerations, 


Asi  " ^vi  " vi-l^  + ^wi  “ wi-l^^  ^ 


sin0 


W • U 

i i-1 


i baj 


COS0 


vi  " Vi-1 


i As , 


The  fundamental  variables  entering  into  the  above  equations  of 
are  defined  as  follows: 

N = internal  axial  force 

Q =*  internal  shear  force 

M = internal  moment 

0 = angle  between  y-axis  and  bar 

As  =•  length  of  bar 

m = discrete  mass 

position  in  the  y direction 


(56) 

(57) 

(58) 

(59) 

(60) 
motion 


v 


w ■ position  in  the  z direction 

■ external  force  per  unit  length  in  the  y direction 

- external  force  per  unit  length  in  the  z direction 

The  double  dot  over  the  variables  v and  w represents  double  differentia- 

2 

tion  with  respect  to  time,  i.e.,  V - d (v)/dt.  Since  each  mass  point 
must  be  in  dynamic  equilibrium  at  all  times,  equations  (55-57)  must 
apply  to  all  points,  i - 1,N,  where  N is  the  total  number  of  masses  in 
the  model. 

The  numerical  procedure  for  solving  equations  (55-57)  is  to 
first  find  all  Q^s  from  equation  (57),  then  find  V . and  from  equa- 
tions (55)  and  (56),  respectively.  Determination  of  N. , M , F , 

yi 

F and  m.  will  be  explained  in  Subsections  4.1.3  to  4.1.8,  and  the 
Zi 

temporal  integration  technique  used  to  extrapolate  quantities  to  the 
next  time  step  is  described  in  Subsection  4.1.10. 

4.1.3  Variable  Cross  Section 

For  beams  of  variable  cross  sections  where  the  width  and 
thickness  of  any  of  the  layers  varies  spanwise  along  the  beam,  a set  of 
cross  sections  is  developed  where  each  corresponds  to  a particular  bar 
in  the  bar-mass  model  (fig.  28).  The  width  and  thickness  is  initially 
specified  at  each  mass  point  and  each  end  of  the  beam.  From  that 
information  uniform  cross  sections  (width  and  thickness  remaining 
constant)  are  generated.  These  uniform  segments  extend  from  one  mass 
point  to  the  next  and  the  idealized  dimensions  of  each  are  found  by 
averaging  the  dimensions  at  either  end. 

Variations  are  assumed  to  be  small  so  that  laterally  applied 
external  pressure  loads  are  assumed  to  act  normal  to  the  beam  as  located 
by  the  mass  points.  In  other  words,  the  local  variations  in  surface 
orientations  are  neglected  in  determining  the  direction  of  the  load 
vector. 
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4.1.4  Idealized  Model  of  the  Cross  Section 


The  creation  of  an  idealized  model  of  the  actual  cross  section 
requires  first  reducing  the  cross  section  to  a set  of  uniform,  rec- 
tangular layers,  as  shown  in  figure  29.  Each  layer,  however,  can  have 
different  physical  properties  as  well  as  different  geometric  properties. 
The  second  step  in  modeling  the  cross  section  is  to  represent  each  layer 
by  a set  of  one  or  more  flanges  in  order  to  represent  the  extensional 
and  bending  stiffness  in  the  layer.  The  more  flanges  used,  the  more 
accurate  the  representation,  but  at  the  expense  of  additional  computer 
time.  Hence,  the  number  of  flanges  and  the  spacing  within  each  layer 
must  be  optimized. 

» 

It  should  be  noted  that  for  a beam  with  nonuniiorm  cross 
section  in  the  spanwise  direction  the  cross  section  at  the  canter  of  the 
beam  is  modeled  for  the  desired  number  of  flanges  and  flange  spacing 
and  this  flange  representation  is  used  for  the  entire  beam.  To  do 
otherwise  would  be  exceedingly  complicated,  resulting  in  little  addi- 
tional accuracy.  This  assumption  does  not  significantly  affect  the 
ability  to  analyze  nonuniform  cross  sections. 

The  number  of  flanges  allotted  to  each  layer  of  material  is 
determined  on  the  basis  of  bending  of  a composite  beam.  The  number  of 
flanges  in  a layer  is  allotted  approximately  in  proportion  to  the  elastic 
modulus,  E,  and  the  width,  b,  for  the  layer,  the  depth,  h,  of  the  layer. 
Because  the  thickness  plays  a more  important  role  in  bending,  the  number 
of  flanges  is  actually  determined  by  h(Eb)^'^.  The  allotment  is  not 
exactly  proportional  to  (h(Eb)^'^)  for  each  layer  since  this  would 
result  in  fractional  numbers  of  flanges.  Actually,  in  the  computer 
program,  the  number  of  flanges  for  each  layer  is  obtained  by  selecting 
an  even  number  n (2,  4,  6,...)  for  the  layer  with  the  largest  value  of 
(h(EB  O'  ^ ) , calculating  a number  n^  for  each  layer  "k"  from 


:(EkV 


0.75 


h(Eb) 


0.75 


(61) 
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and  rounding  n^  Co  Che  nearesC  even  inceger,  or  if  less  chan  1.5,  is 
rounded  Co  1.  Even  integers  are  preferable  because  the  bending  stiff- 
ness of  each  layer  of  material  can  be  represented  more  effectively  with 
a minimum  number  of  flanges.  The  reasons  for  this  selection  are  more 
readily  apparent  from  the  discussions  which  follow  on  the  area  and 
spacing  of  the  flanges. 


The  area  and  spacing  of  the  flanges  in  each  layer  are  evalu- 
ated from  considerations  of  both  the  extensional  and  bending  stiffnesses 
of  the  cross  section.  It  is  desirable,  if  possible,  to  duplicate  the 
bending  and  axial  extensional  stiffness  of  the  entire  cross  section 
through  all  phases  of  stress;  elastic,  elastic  plus  plastic,  and  all 
plastic.  It  is  impossible  to  duplicate  all  possible  elastic  plus 
plastic  stress  conditions  without  the  use  of  an  infinite  number  of 
flanges  in  each  layer  of  material.  It  is  possible,  however,  to  repro- 
duce easily  the  normal  force,  N,  and  the  moment,  M,  resulting  from  the 
following  strain  conditions: 


1.  Any  purely  elastic  strain 

2.  Purely  plastic  axirl  strain 

3.  Purely  plastic  bending  strain 


Satisfying  the  above  strain  conditions  guarantees  that  any  strain  condi- 
tion will  be  reasonably  well  reproduced. 


For  a multilayered  beam,  if  the  elastic  extensional  and  elas- 
tic bending  stiffnesses  of  each  layer  of  material  of  the  cross  section 
are  reproduced,  all  three  of  the  3train  conditions  listed  above  are 
correctly  treated.  Figure  29  illustrates  the  actual  cross  section  for  a 
multilayered  beam  and  the  idealized  flange  representation  of  the  cross 
section  used  in  the  computer  program.  In  the  idealized  representation 
the  flanges  in  a layer  are  of  ‘.qial  area  and  are  distributed  symmetri- 
cally about  the  centerline  of  the  layer.  The  elastic  extensional  and 
elastic  bending  stiffnesses  are  reproduced  by  the  following  area  and 
spacing  of  the  flanges: 
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(62) 


Vk 

“k 


1 


where 


“k 

k 


the  area  of  each  flange  in  layer  "k" 

width  of  layer  "k" 

spacing  between  flanges  in  layer  "k" 

the  height  of  layer  "k" 

the  number  of  flanges  in  layer  "k" 

the  kth  layer 


(63) 


As  pointed  out  earlier,  this  representation  is  based  on  an  even  number 
of  flanges  for  each  layer  of  material.  Actually,  the  uniform  spacing 
between  flanges  could  be  based  on  a symmetrical  distribution  of  flanges 
about  the  centerline  of  the  layer,  with  one  flange  located  on  the 
centerline.  This  representation  is  not  desirable  because  the  flange 
located  at  the  centerline  is  ineffective  in  reproducing  the  elastic  and 
purely  plastic  bending  stiffness  of  the  layer. 

The  one  exception  to  the  spacing  formulas  (63)  is  when  only 
one  flange  is  assigned  to  the  layer,  in  which  case  the  flange  is  located 
at  the  center  of  the  layer  and  ■ 0. 

The  strain  for  pure  bending  is  zero  within  some  layer  of  the 
cross  section.  Therefore,  the  correct  representation  of  the  purely 
plastic  bending  condition  is  obtained  by  dividing  the  layer  in  which  the 
strain  is  zero  into  two  equal  layers.  The  dividing  line  is  the  zero 
strain  axis,  and  both  layers  have  identical  material  properties.  The 
flange  area  and  spacing  are  determined  as  before  except  that  the  ideal- 
ized cross  section  contains  L+l  layers  instead  of  the  L layers  of  the 
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actual  cross  section.  This  approach  is  very  useful  for  a single-layer 
cross  section  because  the  two  layers  are  positioned  symmetrically  about 
the  centroidal  axis.  Representing  the  elastic  bending  stiffness  of  each 
layer  by  flanges  distributed  symmetrically  about  the  centerline  of  each 
layer  reproduces  both  the  elastic  bending  and  purely  plastic  binding 
stiffnesses  of  the  entire  crost;  section. 

The  internal  bending  moment  for  the  cross  section  is  deter- 
mined by  using  the  reference  axis  in  figure  29.  The  inertial  forces  are 
also  based  upon  accelerations  at  the  reference  axis;  therefore,  the 
reference  axis  is  located  at  the  center  of  gravity  of  the  cross  section. 


S(h! + \-i  + r>  Wk 


k-l 


eg 


(64) 


5Z  Wk 


k-l 


where 


eg 


the  location  of  the  centroidal  axis  (reference  axis) 
with  respect  to  the  base  of  the  cross  section 
- the  mass  density  of  layer  "k" 

L - number  of  layers  in  the  cross  section 

and  the  other  variables  are  as  defined  earlier.  As  indicated  in  figure 
29,  the  distance  from  the  center  of  gravity  co  the  rth  flange  located  in 
the  kth  li_yer  is  depicted  by  j;r,  and  represents  the  bending  moment  arm 
for  that  flange.  In  the  case  of  variable  cross  section,  the  parameters 
A^>  d^,  c and  X will  all  depend  on  spanwise  position,  and  henceforin 
the  subscript  "i"  will  be  attached  to  them  to  indicate  that  dependence. 
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A. 1.5 


Strain-Displacement  Relations 


Referring  to  figures  28  and  29  again,  the  finite  difference 
equations  representing  the  strain-displacement  relations  for  the  rth 
flange  in  the  cross  section  at  segment  i are  given  by 


As,  - As, 
i i,o 


N 


As 


i,o 


A9  - A8 , A0,  , - Afl.  . 

i i^o  i-1 i-l,o 


As 


i.o 


As 


i-l,o 


(65) 


'Si  2 


As,  - As. 
i i,o 


As 


i.o 


As.  ~ As... 
i+1 i+l,o 


As 


i+l,o 


where 


, r r (A9.  - A9  ) 

k ^ + rr  ) i k*°- 

2 ui  c’i-t-l' 


As 


i.o 


(66) 


As 


i.o 


2 (Asi  ,o  + Asi+l(o} 


(67) 


A9i  - Vi'0! 


■ sin  (sin8i+^cos0^-cos0i+^sin9:f ) 


(68) 


and 


e “ strain  used  to  compute  the  internal  force,  N. 


= strain  used  to  compute  the  internal  moment,  M. 
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The  change  of  angle  at  mass  i,  A8^,  is  computed 
formula  since  the  sines  and  cosines  are  already 
Likewise,  As^  is  available  from  equation  (58). 
to  the  unstrained  conditions. 


from  the  inverse  sine 
available  (eqs.  (59-60)). 
The  subscript  "o"  refers 


The  two  distinct  formulas  for  computing  and  arise  from 
discretization  - the  fact  that  the  internal  force  is  computed  in  the 
ith  bar;  the  moment  is  computed  at  the  ith  mass.  Once  the  stresses,  a, 
are  calculated  based  on  the  strains,  the  force  and  moment  are  the  result 
of  the  following  summations  over  all  the  flanges: 


(69) 


M. 


i I>;M  k + A^)<;' + ^ 


(70) 


The  following  section  describes  the  stress-strain  relations. 


4.1.6  Stress-Strain  Relations 

The  stress- strain  model,  depicted  in  figures  30  and  31,  con- 
sists of  a set  of  piece-wise  linear  segments  defined  at  coordinates 

((e  ,at),  k = 1,2,. ..n).  This  uniaxial  curve  may  be  different  in  ten- 
k k 

sion  and  compression,  although  the  initial  slope,  is  constrained 

to  be  the  same  in  the  following  mathematic  model. 

In  order  to  incorporate  strain  hardening  and  elastic  unload- 
ing and  reyielding,  a "mechanical  sublayer"  model  is  adopted  (ref.  31). 
This  model  is  based  on  kinematic  hardening  (in  a uniaxial  sense)  which 
takes  the  Bauschinger  effect  into  account  (see  figure  32). 


-93- 


i Wil  in 


ilk 


. YTKHf" 


(a)  Elastic-Perfectly  Plastic  Cycle 


(b)  General  Loading  Cycle 
Figure  32.  Cyclical  Stress-Strain  Load  Paths 
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Each  flange  in  the  program  is  conceptually  broken  down  into 
several  subcomponents,  or  "sublayers".  Each  of  these  sublayers  behaves 
in  an  elastic-perfectly  plastic  manner  and  each  has  the  same  modulus  of 
elasticity,  E,  but  each  has  a distinct  fictitious  yield  stress  associa- 
ted with  it,  The  value  of  is  simply  given  by 


“ Ee, 


(71) 


Once  the  stress  in  each  sublayer  is  found  from  its  own  elastic- 
perfectly  plastic  model  and  the  strain,  e,  the  total  stress  can  be 
determined  by  summing  over  each  sublayer  after  weighting  the  stress 
in  each  sublayer  by  an  appropriate  weighting  factor.  This  weighting 
factor  simply  takes  into  account  the  fictitious  yield  stresses  used 
in  the  model. 


n 


2 ck\<e) 


k-1 


where 


'k 


and 


Wi 


< Wi 


k-1 


k-2,3, . . .n 


(72) 


(73) 


(74) 
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k-n+1 
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In  the  case  of  an  elastic-perfectly  plastic  material,  there  is 
only  one  sublayer.  In  general,  each  distinct  slope  of  the  stress-strain 
model  in  the  entire  tension-compression  regime  will  dictate  a unique 
sublayer,  except  for  zero  slope  segments,  which  are  automatically 
accounted  for.  The  method  makes  two  assumptions:  1)  the  modulus  of 
elasticity  is  the  same  in  both  tension  and  compression,  and  2)  the  break 
point  stresses  increase  for  Increasing  strain. 


For  problems  involving  complex  cyclical  load  paths  where 
points  may  unload,  reyield,  unload,  etc.,  the  method  provides  a very 
convenient  mechanism  to  follow  this  behavior  as  only  elastic-perfectly 
plastic  curves  are  ultimately  involved.  It  is  also  of  note  that  the 


computation  of  the  weighting  factors,  c^'s,  need  only  be  performed  once 


since  they  are  independent  of  time. 


4.1.7  Boundary  Conditions 


DEPROB  has  the  ability  to  analyze  any  combination  of  clamped, 
simply  supported,  and  free  edge  conditions,  including  a symmetric 
"edge",  where  only  half  of  the  beam  or  ring  is  analyzed  when  the  beam 
and  its  external  loading  are  symmetric  about  the  center.  Each  of  the 
edge  conditions  is  discussed  below  and  the  appropriate  equations 
presented. 


Clamped  Edge 


An  ideal  clamped  edge  condition  is  assumed,  where  zero 
deflection  and  zero  slope  are  assumed  at  the  clamp,  as  indicated  in 
figure  33(a). 


According  to  finite-difference  techniques,  the  zero  slope 
constraint  would  best  be  approximated  by  including  at  least  one  mass 
very  near  the  edge,  so  that  the  zero  slope  constraint  is  maintained  in  a 
very  localized  region.  For  large  deflection  problems  dominated  by 


mm 


(a)  Finite  Difference  Model 


(b)  Quasi-Static  Model 


Figure  33.  Clamped  Edge  Models  for  Time  > 0. 
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membrane  forces,  however,  a large  stress  concer c rat ion  occurs  near  the 
clamped  edge  and  results  In  very  large  computational  times  due  to  the 
additional  masses  but,  more  importantly,  also  due  to  the  small  time 
increment  thereby  required  in  the  solution.  The  time  increment,  At,  as 
will  be  seen  in  subsection  4.1.10,  is  proportional  to  the  shortest  link 
in  the  model. 

A new  technique  has  been  applied  which  eliminates  these 
problems,  with  computational  time  on  the  order  of  ordinary  finite- 
difference  calculations,  yet  permitting  a much  higher  degree  of  accuracy. 
This  is  accomplished  by  subdividing  the  segment  (bar)  nearest  the  edge 
into  10  unevenly  spaced  subsegments,  with  the  subsegments  more  closely 
spaced  near  the  edge.  These  subsegments  are  not  used  in  the  finite- 
difference  equations,  however.  Instead  they  are  the  basis  of  a static 
solution,  where,  through  an  iterative  process,  the  forces  and  the  moment 
at  the  end  of  the  first  segment  are  found  at  each  instant  of  time  by 
satisfying  the  end  deflections  and  an  end  slope  condition.  The  use  of  a 
static  solution  is  equivalent  to  ignoring  the  mass  of  the  first  bar, 
which  is  consistent  with  the  discretization  used  throughout. 

Figure  33(b)  shows  the  relative  advantage  of  this  method 
over  the  previous  method  indicated  in  figure  33(a).  The  new  method 
permits  a more  precise  solution  between  the  clamp  and  mass  point  1,  as 
it  is  allowed  to  deform  with  the  rest  of  the  beam. 

When  the  beam  is  initially  modeled,  however,  the  first 
segment  must  be  oriented  perpendicular  to  the  assumed  wall,  as  shown  in 
figure  34(a).  It  is  appropriate  that  the  length  of  that  segment  be  on 
the  order  of  other  segments  in  the  beam.  Either  or  both  of  the  ends  of 
the  beam  can  be  clamped. 


-99- 


(a)  Clamped 


(b)  Simply  Supported 


(N  + i) 

(e)  Free  Ring 


Figure  34.  Edge  Conditions  for  Beam 
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Simply  Supported  Edge 


A no-slip,  knife  edge  (hinge)  condition  is  imposed  by 
requiring  no  deflection  and  no  moment  at  the  edge.  As  shown  in  figure 
34(b),  it  is  appropriate  to  initially  model  the  beam  with  the  segment 
nearest  the  edge  approximately  the  same  length  as  other  segments  in  the 
beam.  Either  or  both  of  the  ends  of  the  beam  can  be  simply  supported. 

Free  Edge 

The  free  edge  is  characterized  by  zero  moment  and  zero 
internal  force  in  the  outer  bar.  As  indicated  in  figure  34(c),  that 
outer  segment  should  be  roughly  half  that  of  other  segments  in  the  beam. 
Either  or  both  of  the  ends  can  be  free. 

In  the  event  that  both  ends  are  considered  free,  the 
structure  can  be  analyzed  if  the  center  of  gravity  remains  stationary 
and  no  rotation  is  introduced.  The  net  translation  and  rotation  of  the 
beam  are  accounted  for  by  calculating  the  accelerations  responsible  for 
this  motion  and  subtracting  them  from  the  left-hand  sides  of  equations 
(55)  and  (56).  The  rigid  body  accelerations  v and  w associated  with 
translation  are  calculated  from 


w 


T 


N 


M 


(75) 


(76) 
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The  rigid  body  rotational  acceleration  is  determined  from 


where 


£ v\ 


i-i 


r^  - (v1  - v)“  + (wi  - w)2 
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(77) 


(78) 


(79) 
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(81) 


(82) 


(83) 


The  coordinates  v,  w represent  the  coordinates  of  the  center  of  gravity 
of  the  system.  Finally,  the  y and  z components  of  acceleration  due  to 
the  angular  acceleration  are,  respectively, 


-9  (w1  - w) 
9 ^vi  ” 


(84) 

(85) 


The  total  rigid  body  motion  is  formed  by  adding  the  contributions  from 
eq.  (75),  (76),  (84)  and  (85). 


Symmetric  Ecge 

When  the  beam  and  its  loading  are  both  symmetric  about 
the  center  of  the  beam,  only  the  first  half  need  be  analyzed.  The  pro- 
gram will  automatically  rotate  the  structure  so  that  the  plane  of  sym- 
metry is  parallel  with  the  vertical  (z)  axis  in  figure  34(d).  A ficti- 
tious mass  (N+l)  is  constrained  to  move  as  a mirror  image  to  mass  N, 
thus  enforcing  zero  slope  and  symmetric  motion  about  the  plane  of 
symmetry.  As  indicated,  the  segment  nearest  the  plane  of  symmetry 
should  be  only  roughly  half  the  length  of  other  segments,  because  its 
mirror  image  is  an  equal  distance  on  the  other  side. 

Free  Ring 

As  indicated  in  figure  34(e),  the  free  ring  must  be 
modeled  symmetrically,  but  is  assumed  to  close  on  itself.  The  program 
then  rotates  the  structure  so  that  the  plane  of  symmetry  is  parallel  to 
the  vertical  (z)  axis.  As  was  the  case  with  a free-free  beam,  rigid 
body  rotation  and  translation  are  subtracted  out  of  the  solution,  leav- 
ing motion  only  relative  to  the  time  dependent  c.g.  of  the  system. 
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4.1.8  Elastic  Support  at  Discrete  Points 

The  flexibility  of  DEPROB,  particularly  the  ability  to  apply 
external  loads  to  discrete  mass  points,  lends  Itself  to  the  analysis  of 
three  special  structural  configurations  of  interest.  They  are  1)  a 
relatively  weak  frame  situated  between  other  relatively  strong  frames 
and  deriving  support  from  the  connecting  longerons  and  stringers,  2)  a 
frame  deriving  support  from  outer  skin  when  internally  loaded  (preblast) 
and  3)  stringers  and  longerons  deriving  support  from  outer  skin  as  in 
(2).  Each  will  be  discussed  below. 

A weak  frame  may  derive  elastic  support  from  a longeron  at  the 
point  at  which  the  two  structures  are  joined.  The  frames  are  assumed  to 
be  evenly  spaced  such  that  the  stiffness  coefficient  associated  with  the 
longeron  is  given  by 

| , 1 (86) 

s i'3 

where 

E'  ■ modulus  of  elasticity  of  longeron 

I'  ■ moment  of  inertia  of  longeron 

V * length  of  longeron  between  strong  frames  or 

bulkheads 

The  frame,  analyzed  by  DEPROB,  will  then  receive  a restoring  force  pro- 
portional to  the  displacement  6 and  applied  to  the  frame  at  that  point. 
This  technique  obviously  excludes  dynamic  coupling  of  the  two  elements 
but  should  be  a reasonable  approximation,  particularly  for  the.  weak 
frame  situation. 

When  a frame  element  is  subjected  to  a net  outward  preload 
pressure  due  to  internal  pressurization  of  the  fuselage,  the  outer  skin 
makes  a very  significant  hoop  stress  contribution  in  resisting  outward 


-104- 


deformation  that  would  not  be  properly  accounted  for  when  modeling  the 
skin  as  part  of  the  frame  cross  section.  In  order  to  correct  for  this, 
the  tensile  stress-strain  curve  is  modified  within  the  program  to  take 
into  account  the  larger  effective  skin  area  whan  the  outer  layer  is  in 
tension. 

For  fuselage  stringers  and  longerons  under  net  outward  preload 
pressure,  a different  approach  is  taken  to  correct  the  same  problem  that 
existed  with  the  frame.  In  this  case,  a cylindrical  section  of  aircraft 
skin  actually  contributes  to  the  hoop  stress  strength  of  the  structure. 
The  restoring  force  at  each  point  along  the  stringer  can  be  shown  to  be 
approximately  proportional  to  the  deflection  of  the  element: 


5 < 0 


-2E0  h sin  6 
s L 


5 > 0 


where 

E * modulus  of  elasticity  of  skin 

h * thickness  of  skin 

9g  * angular  spacing  between  stringers  around  the  fuselage 

4.1.9  Preload  Static  Solution 

The  mechod  chosen  for  solution  of  the  preblast  deformations 
due  to  steady-state  loads  is  an  i.  irative  relaxation  procedure  designed 
to  reduce  the  accelerations  to  zero.  Once  this  equilibrium  state  is 
reached,  the  displacements  are  stored  in  the  program  so  that  a change  of 
blast  range  or  orientation  will  not  necessitate  solving  the  preblast 
problem  again. 


Due  to  the  complexity  of  the  solution  procedure,  the  static 
solution  is  found  for  only  seven  boundary  condition  combinations.  These 
should,  however,  solve  any  practical  aircraft  problem.  The  allowable 
combinations  are 


1st  End 

1)  Clamped 

2)  Clamped 

3)  Clamped 

4)  Clamped 

5)  Simply  supported 

6)  Simply  supported 

7)  Free  ring 


2nd  End 
Clamped 
Free 

Symmetric 
Simply  supported 
Simply  supported 
Symmetric 


4.1.10  Numerical  Analysis 

Integration  of  the  second-order  differential  equations 
(55-56)  is  accomplished  numerically  by  extrapolating  the  2N  displace- 
ments v^,  in  time.  The  integration  used  is  the  central  difference 
formula 

Vi  ■ <ic>2  K + 2-\  - -\-i  <88) 

where 

X represents  one  of  the  2N  displacements 
At  is  the  time  increment 

k denotes  the  time  step,  0,1,2... 


Note  that  step  zero  corresponds  to  static  conditions  just  after  the 
blast  has  arrived  but  before  the  structure  has  responded.  In  order  to 
start  the  procedure,  the  back  point  X ^ must  be  established.  This  is 
done  at  step  zero  by  letting 


(89) 


X - X + 1/2  (At)2  X 
-lo  o 


The  central,  difference  formula  (88)  Is  known  as  an  open 
method.  For  elastic  problems  the  solution  will  be  accurate,  for  all 
practical  purposes,  as  long  as  the  solution  remains  stable.  Too  large  a 
time  increment,  however,  will  trigger  a numerical  instability.  Stability 
can  be  related  to  longitudinal  and  lateral  vibration  frequencies  of  the 
beam  (ref.  30).  These  criteria  depend  on  the  material  properties  and 
the  geometry,  including  the  spacing  between  adjacent  mass  points.  For 
multilayered  beams,  average  material  properties  are  computed;  and  for 
beams  of  variable  cross  section,  each  station  along  the  beam  is  checked 
to  find  the  critical  At.  The  following  equations  represent  the  At 
criteria,  where  the  minimum  value  (when  all  links  i*l,N  are  considered) 
is  appropriate: 


where  n is  the  maximum  number  of  flanges  assigned 
total  thickness  of  the  cross  section,  and 

L 

2 Wk 
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L 
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to  any  layer,  h is  the 


(91) 


The  constant  r has  as  a value  a theoretical  upper  limit  of  1.0,  but  a 
conservative  value  of  0.9  Is  used  to  avoid  marginally  unstable  cases. 

It  should  be  noted,  however,  that  initially  curved  beam  elements  which 
experience  "snap-through"  buckling  may  require  a significantly  smaller 
time  increment  than  that  predicted  by  the  above  formula. 

For  solutions  which  go  inelastic,  and  for  column  buckling 
types  of  problems,  the  solution  may  deteriorate  without  actually  diverg- 
ing. For  these  cases  it  is  recommended  that  particular  attention  be 
paid  to  the  validity  of  the  results.  This  can  easily  be  checked  by 
using  a time  increment  somewhat  smaller  than  that  previously  used  and 
comparing  results. 


4.1.11  DEPROB  Response  Comparisons  with  Experiments  and  with  an 
MIT  Code 


As  partial  verification  of  the  DEPROB  code,  the  dynamic 
responses  of  two  clamped-clamped  uniform  straight  aluminum  beams,  sub- 
jected to  impulsive  loads  on  the  center  portion  of  the  span,  are  pre- 
sented in  table  11  and  figures  35  and  36.  Deflection  time  histories  are 
shown  for  each  beam  as  determined  by  1)  DEPROB  calculations,  2)  exper- 
iment (ref.  29),  and  3)  MIT  beam  code. 


These  particular  beam  tests  were  selected  due  to  the  rela- 
tively good  experimental  data  and  because  of  the  inelastic  nature  of  the 
response.  The  second  beam  is  nearly  identical  to  the  first  except  it  is 
twice  as  thick  and  does  not  exhibit  as  much  strain  hardening  as  the 
first.  In  order  to  make  these  comparisons,  the  DEPROB  program  was 
temporarily  modified  to  incorporate  an  impulsive  lateral  load  over  a 
portion  of  the  beam. 


The  most  important  conclusion  which  can  be  drawn  from  the 
results  is  that  the  two  structural  response  codes  predict  almost  exactly 
the  same  deflection  time  history  (comparison  of  strains  was  not  possible) 
This  is  no  surprise,  since  DEPROB  is  based  to  a large  degree  in  the  MIT 
code;  but,  nevertheless,  enough  changes  have  been  made  so  that  these 
results  are  reassuring.  Secondly,  the  two  experimental  traces  show  no 
significant  differences  from  the  code  predictions,  except  at  late  times, 
thus  lending  even  more  credence  to  the  DEPROB  results. 
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TABLE  11 

GEOMETRICAL  AND  PHYSICAL  PROPERTIES  OF  BEAMS 


DEPROB  Models 

Beam  95 
(ref.  29) 

Beam  111 
(ref.  29) 

Material  Properties 

Aluminum 

Aluminum 

6061-T6 

6061-T6 

E,  psi 

10.81  x 106 

10.7  x 106 

°o»  Psi 

41600 

41200 

Et,  psi 

161000 

61200 

p,  lb/in 

0.098 

0.097 

Length,  in 

10.0 

10.0 

Thickness,  in 

0.124 

0.242 

Width,  in 

1.195 

1.195 

Total  Impulse,  psi-sec 

0.1575 

0.2724 

Initial  Velocity,  ft/sec 

5003 

4480 

Length  of  Charge 
(centered  on  £) 

1.932 

1.988 

Edge  Conditions 

Clamped 

Clamped 

Number  of  Masses 

in  Half  Span  Model 

30 

30 

Number  of  Flanges 

12 

12 

Time  Increment,  sec 

0.75  x 10“6 

0.75  x 10~b 

Measured  Permanent 
Set,  in 

0.581 

0.522 

CLAMPED  BEAM  I0"  X 0.125 


Figure  35.  Dynamic  Response  of  Beam  1/8  Inch  Thick 


CLAMPED  BEAM  10“  X 0. 25 
6061 -T6  ALUMINUM 


Figure  36.  Dynamic  Response  of  Beam  1/4  Inch  Thick 


4.2  DYNAMIC  ELASTIC-PLASTIC  RESPONSE  OF  PANELS  — DEPROP 


The  structural  response  program  DEPROP  was  developed  to  calculate 
the  linear  elastic  and  elastic-plastic  response  of  aircraft  panels  to 
static  and  dynamic  pressure  loads.  Panels  on  aircraft  can  generally  be 
approximated  by  cylindrical  or  flat  panels  with  combination  of  clamped 
and  simply  supported  boundary  conditions  and  can  be  single  or  multi- 
layered with  isotropic  or  orthotropic  material  properties.  The  DEPROP 
program  is  a modification  of  a dynamic  nonlinear  cylindrical  shell 
program  called  DEPICS  (refs.  32-34)  and  contains  the  capabilities  of 
specifying  a static  uniform  pressure  loading  (preload)  and  a transient 
uniform  pressure  loading  on  the  panel. 

The  DEPROP  analysis  is  based  on  the  Novozhilov  nonlinear  strain- 
displacement  relations  for  large  displacement  response  of  thin  panels 
based  on  the  assumption  of  undeformable  normals.  The  program  has 
response  options  for  either  linear  elastic  or  elastic-plastic  material 
behavior.  The  linear  elastic  option  can  be  used  with  multiple  layers  of 
isotropic  or  orthotropic  material.  An  elastically  isotropic  material 
possesses  elastic  properties  which  are  identical  in  all  directions  and 
are,  therefore,  independent  of  the  orientation  of  the  coordinate  axes. 
The  elastically  orthotropic  material  defined  for  these  panels  under 
plane  stress  conditions  possesses  three  orthogonal  planes  of  elastic 
symmetry  that  are  parallel  to  the  geometric  coordinate  axes.  The 
elastic-plastic  option  provides  an. estimate  of  severe  damage  for  a 
single-layered  isotropic  panel  for  a material  with  an  assumed  uniaxial 
bilinear  stress-strain  curve.  The  inelastic  formulation  is  based  on  the 
Mises-Hencky  yield  surface,  a kinematic  hardening  model  and  the  Hencky 
stress-strain  relations  from  the  deformation  theory  of  plasticity  with 
modifications  for  regions  of  elastic  unloading  and  reyielding.  The 
DEPROP  program  calculates  displacement,  strain,  and  stress  time  his- 
tories at  selected  positions  on  the  panel  to  be  used  in  conjunction  with 
various  damage  criteria. 
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4.2.1  Basic  Theory 


The  single-layered  cylindrical  panel  is  considered  to  have  a 
constant  thickness  h,  mean  radius  a,  subtended  angle  9q  and  length  l.  The 
cylindrical  coordinates  (x,  8,  z)  and  the  axial,  tangential,  and  radial 
displacement  components  (u,  v,  w)  are  shown  in  figure  37  on  the  coordi- 
nate surface  which  is  located  at  the  median  surface  of  the  panel.  The 
governing  equations  of  motion  for  the  panel  are  obtained  from  the 
principle  of  virtual  work  for  a dynamic  structural  system  (ref.  35) 
which  is  given  by 


de^dV  - <5T 


<5d  dA 


dt  - 0 


(92) 


where  the  panel  is  undergoing  an  arbitrary  set  of  infinitesimal  virtual 
displacements  6u,  6v,  <5w  that  satisfy  the  geometrical  boundary  conditions 
and  vanish  at  t*t^  and  1 2 » T is  the  kinetic  energy;  0^  are  the  compo- 
nents of  total  stress;  are  the  components  of  total  strain;  F is 
the  surface  force  vector;  d is  the  displacement  vector;  and  integrations 
are  carried  over  volume  V and  deformed  surface  area  a.  It  should  be 
noted  that  this  principle  holds  regardless  of  whether  the  material's 
stress-strain  relations  are  elastic  or  inelastic  and  whether  the  force 
system  is  conservative  or  nonconservative.  If  it  is  assumed  that 
T»T(u,v,w),  then 


3T  ..  3T  ».  3T  - . 
5T  = tv  5u  + — <5v  + tt  ow 
9u  av  3w 


(93) 


where  the  dots  denote  differentiation  with  respect  to  time.  With 
equation  93  and  using  integration  by  parts. 


•5T  dt 


f (-±  ~5u+^l?6v+~r?  owl  dt  (94) 
J 1 dt  3u  dt  cv  dt  ow 
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Cylindrical 


It  is  assumed  that  the  uniform  blast  pressure,  p(t),  acts  on  the 
coordinate  surface  of  the  cylindrical  or  flat  panel.  As  the  panel 
surface  deforms,  the  elemental  pressure  force  vector  remains  normal  to 
the  coordinate  surface  so  that  it  changes  direction  during  deformation. 
The  magnitude  of  this  force  vector  also  changes  as  the  element  surface 
area  of  the  deformed  panel  changes.  It  should  be  noted  that  the  portion 
of  the  pressure  loading  associated  with  the  force  vector's  dependence  on 
the  deformations  represents  a nonconservative  force  system.  Based  on 
the  rectangular  coordinate  system  (X,Y,Z),  the  components  nx»  fly  and 
of  the  inward  unit  normal  surface  vector  and  the  components  dx,  dy 
and  d^  of  the  displacement  vector  d were  defined  in  reference  32  in 
terms  of  u,  v and  w and  their  spatial  derivatives.  Thus,  the  vector  dot 
product  of  the  force  and  virtual  displacement  is  expressed  as 

F • 6d  ■ p(t)(nx<5dx  + HySdy  + nz<5dz)  (95) 

By  neglecting  terms  above  the  second  order  and  recasting  in  terms  of  the 
virtual  displacements  6u,  <5v  and  <5w,  the  virtual  work  done  by  the 
external  forces  is  given  by 


j J F • 5d  dA  - j I p (t)  (N^5u  + Nv  <5v  + 6w)dA 


(96) 


where 


- (w  + w ) 
x x 

O 

-(w  + wg  + v)/a 

O 

1 - (w  + w - vQ)/a  + u 
9 x 

undeformed  surface  area 


The  subscripts  on  the  displacement  components  denote  spatial  derivatives 

O 

and  w denotes  initial  radial  imperfection  in  the  panel. 
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With  equations  94  and  96  and  the  relation 


3e.  3e.  3e 

6e  - ■ ^ 6u  + “ ^ 6v  + -r-^-  6w 

ij  3u  3v  3w 


equation  92  becomes 


3e 


/ (IE  I + fff°  ij  dV  - II  ?Nud*]  5u  + [fe  f 

1 1 TT  A 


J// 


3e 


ii 


ij  3v 


dV 


- //  P»„di]  6v  + jjjj 


3T 


3w 


///■ 


3g  . 


ii 


ij  3w 


dV 


- />><>] 


6w>  dt  = 0 


(9 


The  displacement  components  are  assumed  in  the  following  truncated 


series  form  with  undetermined  time-dependent  coefficients,  u (t), 

mn 


v (t),  w (t): 

mn  mn 


M N 


u(x,e,t) 


z 

m=l  n*»l 


/ umn^x^n  (9) 
7 mn  m n 


M N 


v(x,e,t)  = 


zi 

m=l  n;l 


y (9> 

/ . mn  m n 
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where  $m(x)  and  (0 ) are  functions  that  satisfy  the  geometric  boundary 
conditions  of  the  panels.  The  initial  radial  imperfection  in  the  panel 
is  represented  by 


w(x,0) 


M 

V 

L. 

m*l 


N 


y A ^(x)^(6) 
/ . mn  m n 

n*l 


(99) 


where  A are  prescribed  values  based  on  known  or  assumed  deviations 
mn 

from  the  ideal  shape  of  the  panel.  Based  on  equation  98,  the  following 
relations  are  obtained: 


<5u 


M N 

y ^ Su  4>U4>U»  6v 
mn  m n 

m»l  n*l 


M N M N 

y Sv  <pV4»V  , 6w 
mn  m n ’ 

m*l  n-1  m*l  n“l 


y ) 


mn  m n 


9u 


9u  .u  u 3v 
> 


mn 


m n’  3v 


mn 


J;v  3w  w w 
0_1>  » r- — * « i 
m n 3w  n n 

mn 


( 100) 


Introducing  equation  100  into  equation  97  and  since  6u  , 

mn 

6v||m  and  are  arbitrary,  the  following  3MN  equations  of  motion  are 

obtained: 


3T 

3u 

mn 


d_ 

dt 


+ 


0 


(1013) 


’if  ” 


d_ 

dt 


3T 


oV 

mn 


•///• 


3e 


ij  3 v 


ii 


mn 


dV 


d_  3T 
.dt  3w 

mn 


■///' 


a, 

3e 


ij  3w 


ii 


dV 


mn 


(m-1,2,3. . .M)  (n-l,2,3...S) 


( 101b) 


( 101c) 


; 


i 


f 

p 


'VU  O.V  'UW 

where  the  integrands  of  the  generalized  forces  (QJ  , Q*  , Q^)  are 
given  by 


..  3u  w „ ov 

^ran  P u 3u  ’ P v ov  ’ min  P w 

mn  mn 


3w 

3w 


(102) 


mn 


The  kinetic  energy  of  a single -layered  panel  is  given  as 


l 9o 

T - ^ j j (u2  + v2  + w2)  dxd9  ( 103) 

0 0 


where  o is  the  mass  density  of  the  material  and  the  dots  denote  differ- 
entiation with  respect  to  time.  The  rotary  inertia  contribucions  to  the 
kinetic  energy  have  been  neglected.  Modification  of  the  mass  density 
for  multilayered  pane  s is  introduced  in  subsection  4.2.5.  Further 
development  of  equation  101  depends  upon  the  establishment  of  the 
strain-displacement  relations,  the  stress-strain  relations  and  the 
displacement  component  spatial  functions. 


t- 

; 
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4.2.2  Strain-Displacement  Relations 

The.  strain-displacement  relations  used  in  this  analysis  are 
based  on  the  assumptions:  (1)  strains  are  small  compared  with  unity, 

(2)  the  thickness  of  the  shell  is  small  compared  with  the  radius  and 

(3)  the  Kirchhoff  - Love  hypothesis  that  straight  fibers  which  are 
normal  to  the  undeformed  coordinate  surface  remain  straight  and  normal 
to  the  deformed  coordinate  surface  and  are  not  elongated,  thus  neglec- 
ting transverse  shear  and  normal  strains.  The  basic  formulation  of 
the  following  set  of  nonlinear  strain-displacement  relations  is  attri- 
buted to  Fovozhilov  (ref.  36).  The  total  strain  consists  of  membrane 
and  bending  components  expressed  by  the  fcr^r  ' » e + zk  . The  membrane 
elongation  and  shear  strains  (e  , eQQ,  e ,)  on  the  coordinate  surface 
are  expressed  in  terras  of  the  displacement  components  and  their  spatial 
derivatives: 


. 1 r 2 2 2,  ° 

e * u + -T-  [ w + u + v]+ww 
xx  x 2 x x x xx 


(104a) 


e99  ' « ve  - JXu  + 7T  1(“e  + lv)2 

2a 


+ (v9  - Aw)2  + u02]  +4  Ve 


(104b) 


!*e  ' \ + a u9  + a “x  <“9  + >v)  + a vx(v8 


1 1 ° ° 

+ — uau  + — (w  w,  + w.w  ) 
a 0 x a x 9 9 x 


Aw) 


( 104c) 


Similarly,  the  change  of  curvature  quantities  (ic  , <00.  <x0 ) of  the 
coordinate  surface  which  characterize  the  bending  and  torsional  deforma- 
tions of  the  panel  are  given  by 


(105a) 


K 

XX 


w 

XX 


(1  + vQ/a 


Xw/a  + u ) 
x 


"ee  - ^7  “89  + h ve  + *2  <•“  + v9>  + I \ 

a a a 


+ I3  (w00  + Xve)(v0  ’ Xw) 
a 

, 1 . X . ,2  . X . ,2 

+ ~ ”86ux  + — (ve  - “>  + -jV 

a a a 


(105b) 


+ — w«  (we  + v> 
a 


2 X 2 . 

K m — W + — V + — — W (v  + 3U 

x0  a X0  a x 2 x0  ^ 9 x 

a 


Xw) 


* ^ “x<“9  + V) 
a 


(105c) 


Primarily,  only  those  nonlinear  terms  are  included  in  equation  105  which 
involve  the  radial  displacement  and  its  derivatives.  The  subscripts  on 
the  displacement  components  in  equation  104  and  105  denote  partial 
spatial  derivatives.  The  end  terms  in  equation  104  are  included  to 
account  for  the  initial  radial  imperfection  of  the  panel  as  indicated  by 
Donnell's  representation  in  reference  37.  The  parameter  X is  introduced 
in  the  strain-displacement  relations  so  that  they  apply  to  both  curved 
and  flat  panels.  Thus,  X = 1 for  curved  panels.  For  flat  panels, 

X * 0,  a = 1,  0 is  replaced  by  y,  and  9q  is  replaced  by  b,  the  width 
of  the  flat  panel. 
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4.2.3  Constitutive  Relations 


In  DEPROP,  the  behavior  of  the  panel  material  is  treated  as 
elastic-plastic  for  Isotropic  single-layered  panels  and  elastic  for 
isotropic  and  orthotropic  multilayered  panels.  The  elastic-plastic 
analysis  for  the  single-layered  panel  has  been  established  as  the  basic 
formulation  in  the  DEPROP  program.  The  elastic  multilayered  analysis  is 
established  as  an  alternate  option  based  on  appropriate  modifications  of 
the  elastic-plastic  formulation.  In  the  DEPROP  analysis  the  solution 
involves  total  strains  and  stresses;  therefore,  for  response  in  the 
inelastic  region,  it  is  convenient  to  use  the  deformation  theory  of 
plasticity  instead  of  flow  theory  which  involves  incremental  strains  and 
stresses.  Plastic  deformation  theory  is  based  on  an  averaging  process 
that  permits  a total  strain  solution  dependent  upon  only  the  final 
stress  state  at  the  end  of  a loading  path.  In  general,  deformation 
theory  is  an  approximation  of  the  more  rigorous  flow  (incremental) 
theory  but  is  equivalent  to  flow  theory  for  an  elastic-plastic  material 
when  the  stress  loading  is  proportional,  that  is,  the  ratio  of  principal 
stresses  remain  constant  during  the  loading  process.  However,  since  the 
dynamic  response  solution  is  solved  incrementally  in  time  by  numerical 
methods  in  DEPROP,  the  strain  increments  are  small  and  the  stress  state 
is  fairly  constant  in  the  plastic  region  over  each  time  step  for  which 
the  equations  of  motion  are  solved.  Thus,  the  plastic  deformation 
theory  provides  a much  more  accurate  solution  when  the  averaging  process 
takes  place  separately  over  each  small  time  increment  as  the  response 
solution  is  obtained  by  a step-by-step  timewise  procedure. 

f1 

In  deformation  theory  the  total  strain  is  a function  of  the 
state  of  stress  and  consists  of  a recoverable  elastic  component  and  a 
nonrecoverable  plastic  component.  It  is  assumed  that  the  material  is 
incompressible,  that  is,  no  permanent  change  in  volume,  due  to  the 
plastic  strain.  Thus,  the  total  plastic  strain  is  equal  to  the  deviator 
plastic  strain.  Furthermore,  it  is  assumed  that  the  material's  uniaxial 
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stress-strain  curve  is  modeled  by  the  bilinear  representation  shown  in 
figure  38  in  which  the  strain  hardening  is  defined  by  slope  Et<  This 
stress-strain  representation  is  interpreted  for  the  biaxial  state  of 
stress  through  the  use  of  the  effective  stress  (o)-effective  strain  (e) 
concept  in  which  the  secant  modulus  (Eg)  indicated  in  figure  38  is 
defined  by 


E* 

s 


a + E„  (e  - e ) 
o t o_ 

G 


(106) 


where  a , gq  are  the  yield  stress  and  strain,  respectively,  from  the 
material's  uniaxial  bilinear  representation.  Thus,  the  effective 
stress,  effective  strain  and  secant  modulus  quantities  are  used  to 
relate  the  biaxial  stress-strain  condition  to  the  assumed  uniaxial 

bilinear  stress-strain  representation  for  the  isotropic  material.  The 

— — *\, 

effective  stress  and  strain,  expressed  as  u * f(a^)  and  e * gU^)* 
are  functions  of  the  total  stress  and  strain  components,  respectively, 
and  are  more  conveniently  introduced  in  explicit  form  later  in  the 
development . 

The  Hencky  stress-strain  relations  for  deformation  theory 
(ref.  38)  are  used  in  the  plastic  region  and  are  given  in  the  following 
form: 


ij 


I t(l+n)«lj-w>1[k51jl 


+ 1 <r  - 

s 


(107) 


where  E is  the  modulus  of  elasticity,  v is  Poisson's  ratio  and  5 is  the 
Kronecker  delta.  The  first  portion  of  equation  107  represents  the 
elastic  component  of  strain  while  the  second  portion  represents  the 
plastic  component  of  strain  in  terms  of  the  deviator  stress.  For  use  in 
equation  101,  the  stress-strain  relations  in  equation  107  are  inverted 

into  the  form  a. 


ij 

and  are  given  by 


f (e..)  for  the  case  of  plane  stress  (a  =o3  =o  =0) , 
ij  r zz  9z  xz 


i 

4 
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+ vAk5ij] 


(i»J  »k“l,2) 


(108) 


where 


1 s .1  . 

va  ‘ 2 ' r <2  - u) 


e12  " 2 ex6’  £11  “ £xx’  e22 


It  should  be  noted  that  the  forms  of  the  stress-strain  relations  in  the 


elastic  (£  »E)  and  plastic  regions  are  identical, 
s 


The  yield  criterion,  in  conjunction  with  a hardening  rule,  and 
the  stress-strain  relation  for  unloading  and  reyielding  by  which  the 
past  strain  history  is  preserved  are  to  be  established.  The  initiation 
of  yielding  for  a biaxial  state  of  stress  is  based  on  the  Mises-Hencky 
yield  criterion  (ref.  38)  given  as 


t0ll  + °22  " all  a22  + 3°12  J 


(109) 


where  a is  the  equivalent  or  effective  stress  and 


°11  " axx’  °22  " a99  * °12 


This  yield  criterion  states  that  plastic  flow  will  occur  when  the  equiv- 
alent stress  a reaches  a value  equal  to  the  uniaxial  yield  stress  in 


tension  oq.  A kinematic  hardening  model  is  employed  in  conjunction  with 


the  Mises-Hencky  yield  surface  which  accounts  for  the  Bauschinger  effect 
when  reyielding  occurs  due  to  the  strain  reversals  during  unloading. 

The  Bauschinger  effect  for  a strain  hardening  material  is  described  by 
the  yielding  behavior  of  a material  at  a reduced  yield  stress  when 
reloaded  in  the  opposite  direction  from  that  of  the  initial  yielding. 


.-A 


The  kinematic  hardening  models  discussed  in  reference  39  assume  chat 
during  plastic  deformation  the  yield  surface  translates  as  a rigid  body  in 
stress  space  with  the  size,  shape  and  orientation  of  the  elliptical  yield 
surface  being  invariant.  The  kinematic  hardening  model  to  be  used  in 
this  analysis  is  illustrated  in  figure  39  for  the  Mises-Hencky  yield 
surface  in  the  plane  of  the  principal  stresses  and  o^.  Corresponding 
to  the  initial  yielding  position  (i)  and  the  unloading  position  (f) 
indicated  in  figure  38,  the  rigid  translation  of  the  yield  surface  for  a 
shift  of  the  stress  state  from  position  (i)  to  position  (f)  is  shown  in’ 
figure  39.  The  change  in  total  stress  components  from  position  (i)  to 
position  (f,  are  defined  by  and,  similarly,  the  corresponding  change 
in  the  total  strain  components  are  defined  by  ^ , so  that 

'v  r 'vr-l  , r-1  r-1 

“ij  " aij  + aij(f)  " aij(i) 


(110) 


■vr-1  'vr-l  ^r-1 
Bij  + eij(f)  “ £ij(i) 


where 


r ■ the  number  of  elastic  unloadings  from  yielded  conditions 
(r-1, 2...) 


a°  - T - 0 

^ x ± **  J a 


ij  "ij 

(i)  indicates  initiation  of  yielding  or  reyielding 
(f)  indicates  final  position  prior  to  unloading 


The  yield  criterion  for  the  translated  yield  surface  is  based  on  the 
effective  stress  given  as 


■>r  . 2 


a - [(au  «u)  - (°11-ai1)(°22"a22)  + (a22“^22)2  +3(a12-a"l2)2] 

(111) 

Furthermore,  it  is  advantageous  in  this  analysis  to  relocate  the  origin 
on  the  e axis  after  each  unloading  such  that  the  extended  elastic 
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unloading  curve  passes  through  the  zero  position.  This  is  accomplished 
by  defining  the  effective  strain  as  follows: 


t - ) 


[^2  [<1'%'fV5)  Kl  - ^1>2  + <=22^2>2) 


(112) 


6,,-2rr„) 


a+v  ,2  ” 


r 3)  1/2 


Thus,  the  elastic-plastic  behavior  of  the  material  for  sub- 
sequent yieldings  after  an  unloading  has  occurred  is  always  based  on  the 
same  a versus  e curve  which  originates  at  position  (0,0).  This  approach 


requires  that  the  stress-strain  relations  be  modified  by  the  and 
?^4  quantities  for  unloading  and  reyielding  conditions  to  account  for 


the  past  stress-strain  history.  The  general  fora  of  the  stress-strain 
relations  for  the  elastic,  elastic-plastic,  elastic  unloading  and 
plastic  reyielding  regions  are  identical,  so  that  the  general  stress- 
strain  relations  based  on  the  form  of  equation  (108)  is  given  by 


^r 

“u  * 


+ (u3> 


(i,j,k  - 1,2) 


where  for  the  following  regions  of  response, 


a)  initial  elastic  loading 


E =E,  £.r  = 3.r  - 0 
s ij  ij 


b)  initial  plastic  loading 


'v#  r r 

E *E  , a . . * 3 = 0 

s s’  ij  ij 


c)  qth  elastic  unloading 


d)  qth  reyielding 


_ P ^ r 
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Thus,  there  are  four  basic  regions  of  response  for  which  the  stress- 
strain  relations  have  been  established  by  equation  113.  For  an  elastic- 

'll 

perfectly  plastic  material,  E^  ■ 0 and  a are  set  equal  to  zero  in 
equations  111  and  113.  It  should  be  noted  that  for  a strain  hardening 
material,  a stress  path  which  may  move  along  the  yield  surface  (neutral 
loading)  would  not  be  properly  represented  in  the  analysis,  since,  upon 
unloading,  the  yield  surface  would  be  rigidly  translated. 

For  elastic,  isotropic  or  orthotropic  multilayered  panels,  the 
stress-strain  relation  formulation  follows  the  approach  presented  in 
reference  40.  In  orthotropic  layers,  the  geometric  cylindrical  coordi- 
nate axes  and  principal  orthotropic  direction  are  assumed  parallel.  The 
multilayered  cross  section  for  the  panel  is  shown  in  figure  40  with  the 
nomenclature  used  in  the  following  formulation.  The  position  of  the 
coordinate  surface  relative  to  the  inner  surface  of  the  panel  is  defined 
by  the  distance  H.  The  membrane  and  bending  stress  resultants  for  the 
multilayered  panel  are  given  by 


°xx  * Cllexx  + C12£0e  + F11<XX  + F12Kee 
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ae«  " c22£ee  + ci2Exx  + F22icee  + Fi:K 


XX 
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axx  “ DllKxx  + D12<9d  + F11EXX  + F12e90 


J90  = °22<90  + °12<xx  + F22£00  + F12Exx 


ax6  = D33<x6  + F33£x9 


(114a) 
(114b) 
(114c) 
(I14d) 
( 1 14e ) 
( liar ) 
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by 


The  elastic  stiffness  constants  C..,  F . . and  D..  are  defined 
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where 


ir  k 

E , E„  are  the  moduli  of  elasticity  in  the  x and  9 directions, 
x 9 


respectively,  of  the  kth  layer 


k k 

v , V-  are  Poisson  s ratios  in  the  x and  9 directions,  respec- 


x9 


tively,  of  the  kth  layer 

is  the  shear  modulus  of  the  kth  layer 

is  the  distance  from  the  inner  shell  surface  to  the  outer 
surface  of  the  kth  layer 

is  the  total  number  of  layers 


For  an  isotropic  material  E = E,  - E,  v - v,  = v and  G , = 7 . 

x 9 x 9 x0  2(l+v) 
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Ic  has  been  found  chat  Che  optimal  position  of  the  coordinate 
surface  for  the  most  efficient  modal  convergence  is  at  the  neutral  axis 
of  the  cross  section.  When  the  coordinate  surface  is  located  at  the 
neutral  axis,  the  interaction  stiffnesses  vanish.  These  interaction 
stiffnesses  reflect  the  influence  of  the  change  in  curvature  on  the 
membrane  stress  resultants  and  the  membrane  strains  on  the  bending 
stress  resultants.  As  discussed  in  reference  40,  for  the  general  case 
of  an  antisymmetrical  orthotropic  multilayered  cross  section,  a neutral 
axis,  which  is  defined  when  all  F ■ 0,  does  not  exist  except  for 
special  combinations  of  elastic  characteristics  of  the  various  layers. 
For  the  general  case,  the  position  of  the  coordinate  axis,  defined  by  H, 
is  established  for  this  analysis  by  setting  F^  ■ F00  * * F^  * 0 to 

obtain  the  values 


Z -u  (h*  - 'Vi) 

H - ^ (115) 

ij  n 

2 E Bij (h*  - iw 

k»l 

and  then  H is  determined  by  averaging  these  values  as  follows: 

B - ? (SU  + *22  + *12  + a33>  (li4) 

It  should  be  noted  that  for  cases  where  the  neutral  axis  does  exist,  the 
coordinate  surface  is  located  at  this  position  through  the  above  proce- 
dure. When  the  center  of  mass  of  the  cross  section  does  not  coincide 
with  the  neutral  axis,  a slight  discrepancy  in  the  inplane  inertia  would 
be  introduced  since  the  rotary  inertia  is  not  included  in  this  analysis. 
Experimental  results  indicate  that  the  rotary  inertia  affected  the  response 
quantities  by  about  only  1Z  for  a shell  undergoing  large  displacement  response. 


i 


a 

1 


i 
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<1 

i 

{ 

\ 
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4.2.4  Displacement  Component  Functions 

In  equation  98  the  displacement  components  are  expressed  in 

series  form  as  a product  of  time-dependent  coefficients  and  independent 

spatial  functions  <p  (x)  and  t (9).  These  spatial  functions  are  selec- 
m n 

ted  so  as  to  satisfy  the  geometric  boundary  conditions  of  the  panels. 
The  boundaries  of  the  panel  are  assumed  to  be  either  clamped  or  simply 
supported  and  spatial  functions  are  defined  to  cover  all  combinations 
of  these  boundary  conditions  for  the  four  edges  of  a panel  defined  by 
x * 0,  i and  9 - 0,  9 . On  clamped  edges  the  boundary  conditions 


w ■ v 


u 


5w  _ _3w 

3x  9£ 


0 


are  to  be  satisfied  while  on  simply  supported  edges  the  boundary  condi- 
tions 


2 2 


are  to  be  satisfied.  Since  the  panels  are  uniformly  loaded,  the  assumed 

displacement  functions  will  adhere  to  this  symmetry.  The  nondimens ional 

TTX  ^0 

variables  y * — and  8 * - — are  introduced  for  use  in  this  analysis. 

The  spatial  functions  for  ?he  u and  v displacements  are  assumed  to  be 
the  same  whether  thn  edges  are  clamped  or  simply  supported  and  are 
given  by 


4>U(y)  * sin  2my 

?U(3)  = sin  (2n-l)S  (117) 

n 

4>V(y)  = sin  (2m-l)y 
m 

$V(8)  = sin  2n8 
n 
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The  boundary  combinations  for  the  y and  3 directions  are  based 
on  opposite  edges  being  both  clamped,  both  simply  supported  or  one 
clamped  and  one  simply  supported.  The  w-displacement  functions  for  the 
y and  8 directions  are  based  on  the  natural  vibratory  mode  shapes  of  a 
uniform  beam  with  corresponding  end  boundary  conditions.  These  spatial 
functions  in  the  y and  8 directions  are  given  as  follows  for  the  three 
boundary  combinations: 

For  clamped/clamped  or  clamped/simply  supported 


w 

* 

■ 

1 

s 

*<  \ 

cos 

A Y 

m 

m 

cosh 

TT 

TT 

w 

cosh 

X s 

n 

X 3 

n 

i a 

n 

TT 

cos 

TT 

X Y 

sinh sin 


a«  (' 

(/\  p 

sinh 

TT 


X_Y 


IT 

X 3 


- sin 


f) 


(118) 


where 


X or  X are  the  odd  roots  of  cos  X,  cosh  X.  ■ 1 for  the  clamped/ 
m n i i r 

clamped  boundary  condition 

X or  X are  the  odd  roots  of  tan  X,  = tanh  X.  for  the  clamped/ 
m n li 

simply  supported  boundary  condition 


a . 


cosh  X.  - cos  X. 

i i 


i sinh  X,  - sin  X , 
l i 


n or  m 


For  simply  supported/simply  supported 
w 

a 3 sin  (2m-L)y 
m 

;W  = sin  (2n-l)3 


(1191 


It  should  be  noted  that  the  functions  given  in  equations  118  and  119  are 
orthogonal . 
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For  the  general  curved  panel,  there  are  nine  combinations  of 
boundary  conditions  provided  by  the  DEPROP  program,  namely, 


y-direction 


6-direction 


The  first  four  combinations  (1-4)  have  symmetry  in  both  directions,  the 
next  four  combinations  (5-8)  have  nonsymmetry  in  one  direction  and 
combination  9 has  nonsymmetry  in  both  directions. 


4.2.5  Governing  Equations  of  Motion 

With  the  strain-displacement  relations  (equations  104  and 
105),  the  stress-strain  relations  (equations  113  and  114)  and  the  dis- 
placement component  functions  (equations  117-119)  defined,  the  governing 
equations  of  motion  (equation  101)  for  elastic-plastic  deformations  are 

developed  further  by  performing  the  indicated  spatial  integrations.  For 

W V u 

convenience,  the  dimensionless  quantities  W » — , V = — , U = — , 

o ° - , 3 d 3. 

tT  W 7TX  _ t i . Tr  n a , v 

W , Y = — , 3 = — , L = — , J = — , R = ^ and  K = ca  are 

introduced  into  the  formulations.  With  this  notation  and  the  spatial 
integration  for  the  kinetic  energ’  i equation  103  performed  analyti- 
cally, the  governing  equations  of  motion  in  the  radial  displacement 
direction  (equation  101c)  are  given  by 
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IT  TT  h/2 


\Vt2“mn  + f-  / / I 


3e 


-h/2 


ij  aw 


ii 


dydBdz 


mn 


-III' 


dydB 


(m-l,2,...M) 
(n-1,2, . . .N) 


( 120) 


where  for  the  w-equations 

k , k - /T  for  C-C  or  C-S  opposite  boundaries 
1 P 

= 1//T for  S-S  opposite  boundaries, 
for  the  u and  v-equations 

k - k - 1//T 

Y P 


and 


'v  9 0 I iu 

% - 2LZRp(i  - XW  - XW  + JVg  + ± Uy) 


mn 


(121) 


Qv  = -2L2Rp(JW3  + JW3  + XV)  -rf- 


mn 


Q = -2LRp  (W  + W ) ... 

u 7 7 aU 


au 


am 


Although  equation  120  is  given  in  terras  W , the  equations  of  motion  for 

the  tangential  and  axial  displacement  directions  can  be  obtained  by 

replacing  W with  V and  U , respectively,  and  using  the  appropriate 
mn  mn  mn 

k , kg  and  Q expressions. 

The  remaining  spatial  integrations  in  equation  120  are  to  be 
accomplished  numerically  thus  providing  a mechanism  for  discretization 
through  the  spatial  points  selected  to  compute  the  representaj:!^^ 
elastic-plastic  behavior  throughout  the  panel.  Thus,  a sufficient 
number  of  spatial  points  must  be  specified  to  obtain  a satisfactory 
deformation  response  solution.  For  integration  through  the  thickness  of 
the  panel  in  the  z direction,  it  is  convenient  to  separate  the  integrand 
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Into  parts  which  either  are  or  are  not  explicitly  dependent  on  the  z 
variable,  that  is,  involving  membrane  strains  and  bending  strains.  The 
total  strain  quantities  given  in  equation  120  for  an  arbitrary 
position  in  the  panel  consist  of  the  membrane  and  bending  components 
given  by 


*11  ’ £u  + 


(122) 


Therefore,  the  integrand  can  be  given  by  fm  + zf*3  where 


..m 

f - a 


3e 


xx 


xx  3W 


+ 0 


mm 


3e 


60 


06  3W 


mn 


3e 


+ a 


x6 


x6  3W 


mn 


(122a) 


, 3K 

cb  xx 

f - a 


xx  3W 


mn 


3K 


30 


00  3W 


mn 


3K 


x6 


x0  3W 


mn 


(123b) 


and  the  total  stress  components  are  obtained  from  equations  113 
and  110  in  which  i,j=l  denotes  x and  i,j=2  denotes  0.  The 
Legendre-Gauss  quadrature  formula  (ref.  41)  was  chosen  for  the 
numerical  integration  in  the  z direction  where  L is  the  number  of  points 
selected  through  the  thickness  of  the  panel.  In  the  y and  S directions 
it  is  convenient  to  have  even  spacing  and  it  is  advantageous  to  have 
spatial  points  on  the  clamped  edges  and  at  the  center  of  the  panel. 
Simpson's  quadrature  formula  (ref.  41)  satisfies  these  desirable  features 
and  therefore  was  selected  over  various  Gaussian  quadrature  formulas. 

The  number  of  spatial  points  selected  in  the  y and  8 directions  are 
given  by  M and  N,  respectively,  where  M and  N must  be  odd  numbers.  By 
performing  the  indicated  numerical  integrations,  equation  120  is  recast 
into  the  form: 
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k k0pi2W__  + 


y 6 mn 


9 (M-l) (N— 1 ) 


I £ Vm 


j-l  k-l 


L2  £ Ht  [^(Yj.V 


i-1 


+ ~2R  5ifi  (vM  (Yr  v}-° 


(m-l,  2.. .M) 
(n-1,2. . .N) 


where 


5 are  the  zeros  of  the  Legendre  polynomial  P-  (£) 
X.  L 


ll  r 

2i  * 2 


H 


2(i-eJ) 


1 d>»2  (pj^tEpi2 


H . - 4 (j  ,k  - even) 
j or  k J 


2 (j,k  - odd,  except  for  j-l,  M or  k - 1,N) 
1 (j-l,  M and  k-l,N) 


^ '(£) ' 


'N  “1  / 


(124) 


When  symmetry  is  present  in  both  the  y and  3 direction,  for  example, 
only  one  quarter  of  the  panel  need  be  considered  and  the  spatial  inte- 
gration takes  the  form 


-137- 


it/2  tt  / 2 
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where 


c 


J 


- 1 


j ■ 1, 


2. 


M-l 

2 


1 

2 


d 


k 


- 1 


k - 1,  2, 


N-l 

2 


1 

2 


k - 


N+l 
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For  the  purely  elastic  solution  of  a panel,  the  integrations  through  the 
thickness  can  be  obtained  analytically,  and  results  in  the  following 
simplification  in  equation  124: 


where  replaces  a In  fm  (equation  123a),  replaces  in 

f*3  (equation  123b)  and 


ij  ’ x_y2  L(1'u)£ij  + V£kk5ij] 


(127a) 


(i,j,k  - 1,2) 


(127b) 


The  quantities  e^,  eQQ  and  exQ  a » given  by  the  normalized 
versions  of  equations  104  and 
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The  quantities  K , K and  K are  given  bv  the  normalized  versions  of 

AA  v XU 

equations  105  and 
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mn 


For  the  elastic  response  solution  of  multilayered  panels,  the 

same  formulation  is  used  with  several  modifications.  The  stiffness 

constants  in  equation  114,  C.  , F and  D,  , are  considered  to  have  been 
2 3 *3 

divided  by  a,  a and  a , respectively,  and  R * a/h^.  equation  124, 

P is  replaced  by  p given  by 


k 2 


\ - w 


k-l 


and 


Hif” 


2Rf‘ 


i=*l 


4R2£b 


i-1 


The  am  and  a*3  quantities  given  in  equation  114  replace  the  appropriate 
total  stress  quantities  given  in  equation  123.  With  these  modifications 
in  DEPROP,  elastic  response  solutions  can  be  obtained  for  multilayered 
panels  of  isotropic  or  orthotropic  material  layers. 
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4.2.6  Numerical  Analysis 

The  second-order  differencial  equations  given  by  equation  124 

and  corresponding  equations  for  V and  U are  to  be  solved  numerically 

nm  mn 

in  time.  The  integration  method  used  to  obtain  an  approximate  timewise 
step-by-step  solution  is  based  on  the  central  difference  formula  given 
by 

Vi  ■ K <4t)2  + 2\  - \-i  <i3o> 

where 

X represents  the  normalized  undetermined  time-dependent  displace- 
ment coefficients,  W , V and  U 

rap  mn  rap 

At  is  the  time  increment 
k denotes  the  time  step 

This  central-difference  method  replaces  the  higher  order  integration 
method  used  in  the  DEPROP  program  in  reference  1.  The  central-difference 
method  permits  a time  increment  about  50%  greater  than  the  previous 
higher  order  integration  method  without  any  significant  change  (less 
than  1%)  in  the  response  solution  for  the  panel  problems  considered. 

The  starting  procedure  used  for  this  method  is  the  same  as  that  pre- 
sented in  subsection  4.1.10. 

In  solving  the  set  of  simultaneous  second-order  differential 
equations,  spatial  integrations  must  be  performed  in  the  y and  3 direc- 
tions and  in  the  z direction  for  the  elastic-plastic  solution  during  the 
stepwise  time  integration.  The  required  integrations  are  performed 
numerically  during  each  time  step  using  the  values  of  the  displacement 

coefficients  W , V and  U for  the  particular  time  step  to  compute 
mn  mn  mn 

the  displacements  and  their  derivatives,  the  strain  quantities  and  the 
stress  quantities  used  in  equation  124. 
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Several  situations  arise  in  the  Implementation  of  the  biaxial 
elastic-plastic  theory  in  DEPROP  which  require  special  numerical  treat- 
ment. These  situations  are  associated  with  the  overshoot  during  the 
time  increment  in  which  yielding  occurs,  the  criteria  for  determining 
elastic  unloading  and  restrictions  if  unloading  is  followed  by  immediate 
reyielding,  and  the  consistent  determination  of  e,  Es>  and  during 
each  time  step.  The  special  numerical  schemes  used  to  treat  these  three 
situations  are  described  briefly  in  the  following  three  paragraphs, 
respectively. 


Whenever  a point  in  the  panel  yields  or  reyields  during  a time 
increment  (o  > 0Q) . Che  stresses  cannot,  in  general,  be  computed  on  a 
purely  elastic  basis.  The  computation  of  stresses  should  follow  the 
bilinear  stress-strain  curve;  but  this  is  very  difficult  to  effect  since 


o is  not  a linear  function  of  the  o 


iJ 


Instead,  an  iterative  scheme 


is  employed  to  adjust  the  o^'s  proportionately,  so  that  o - oq.  By 
elastic  relations,  the  associated  strains  (o^)  are  then  determined 
for  the  later  computation  of  (equation  1J0).  The  process  for 
correcting  for  overshoot  when  yielding  occurs  between  times  t^_^  and  t^ 
is  illustrated  in  figure  41.  The  values  of  e and  o are  shown  at  time 
t^  the  values  indicated  at  t^  represent  hypothetical  uncorrected 
values.  By  linearly  adjusting  the  stresses,  the  point  (e^,  c ) is 
reached.  It  is  noted  that  the  actual  point  at  time  t^  should  be  (e  , 
oA)  instead  of  (e^,  oq),  but  the  error  in  stress  is  relatively  small 
since  Et/E  <<  1.  The  error  introduced  is  proportional  to  the  size  of 
the  integration  time  step  used. 


For  a point  in  the  panel  which  is  yielding  at  time  t,  . , elastic 

unloading  is  detected  when,  in  proceeding  to  the  next  time  r,  , the  equiva 

lent  strain  decreases,  i.e.,  < c.  When  this  occurs,  £,  and  j,  are 

it  k-1  k k 

recomputed  using  the  elastic  unloading  version  of  equation  113.  Fur- 
thermore, unless  o^  is  less  than  , it  is  assumed  that  the  point  did 
not  unload.  This  possible  inconsistency  is  partially  numerical  in 
nature,  and  is  partially  due  to  the  nonlinearity  or  the  equations 


Figure  41.  Correction  for  Overshoot  at  Yielding 

111  or  112.  However,  only  rarely  will  a point  pass  the  strain  criterion 
for  unloading  but  fail  the  stress  criterion.  Due  to  numerical  discrep- 
ancies, it  is  possible  for  the  computed  in  the  unloading  region  to  be 

tC 

greater  than  0q  at  time  t,^.  This  inconsistency  results  in  reyielding 
without  any  overshoot  correction  being  made.  In  fact,  this  event  repre- 
sents a numerical  error  and  is  usually  associated  with  the  initial 
stages  of  numerical  instability  of  the  solution.  Consequently,  if  this 
event  frequently  occurs,  the  run  is  automatically  terminated  and  a 
smaller  time  increment  must  be  selected. 

In  the  temporal  integration  sequence,  the  displacement  coefficients 
are  computed  for  the  end  of  the  next  time  step  at  t,  through  the 

tv ' X 

central-difference  formula  (equation  130)  given  the  past  aisplacement 

coefficients  at  t,  and  t,  , and  the  acceleration  at  t,  . These  extrao- 
K.  K-l  K. 

olated  displacement  coefficients  are  then  used  to  compute  . . at  t. 

l ] <t+l 

Then,  for  points  in  the  plastic  region,  the  quantities  , E and  are 

S 3 
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evaluated  in  order  to  compute  a at  t^.  However,  e is  dependent  on 
Vg  (equation  112),  Eg  is  dependent  on  e (equation  106),  is  dependent 
on  E (equation  108),  and  a is  dependent  on  E and  v (equation  113). 
Thus,  a simple  iteration  scheme  is  used  to  simultaneously  solve  for  the 

three  parameters  e , E and  v , so  that  a consistent  determination  of 

s s 

°ij  c,n  ke  ma<*e  at  Ck+1‘  ^en»  wich  the  displacements,  their  derivatives, 
£ij  anc*  °ij  » determined  at  t^+^,  the  accelerations  are  computed  at 
t^+^  through  equation  124  and  the  whole  process  is  then  repeated  for  the 
next  time  step. 


A method  has  been  established  for  estimating  time  increments  it  for 
the  temporal  integration  that  would  result  in  stable  solutions  for  the 
majority  of  panel  cases.  The  proper  time  increment  is  a function  of 
geometric  and  physical  properties  of  flat  or  curved  panels  and,  for  the 
DEPROP  formulation,  also  a function  of  the  number  of  modes  used  and  the 
spacing  between  spatial  integration  points.  The  method  for  estimating 
time  increments  is  based  on  formulas  for  the  higher  vibratory  frequen- 
cies of  linear  elastic  panels  which  incorporate  the  aforementioned 
parameters.  The  basic  frequency  formulas  for  single-layered  flat  and 
curved  panels  were  obtained  from  references  42  and  43  and  modified  for 
multilayered  panels.  The  time  increment  is  estimated  by  the  product  of 
the  reciprocal  of  the  frequency  and  an  arbitrary  adjustment  factor.  For 
both  threshold  of  permanent  and  catastrophic  damage  levels,  the  panel 
response  is  nonlinear  so  that  the  arbitrary  adjustment  factors  are 
determined  by  back-figuring  from  the  time  increments  found  to  give 
stable  solutions  for  various  representative  panels.  In  some  panel  cases 
considered  where  the  spacing  between  integration  points  was  critical 
for  numerical  stability,  it  was  found  that  the  it  formula  used  in  refer- 
ence 44  for  finite  difference  solutions  is  applicable.  This  formula  is 
based  on  the  time  for  an  inplane  elastic  wave  to  propagate  between  mesh 
points . 


For  flat  panels  the  governing  time  increment  used  as  the  initial 
estimate  for  it  in  DEPROP  solutions  is  the  smaller  it  obtained  from  the 


following  two  formulas.  The  first  formula  is  associated  with  the 
highest  solution  frequency  of  a flat  panel  with  inplane  stresses  which 
are  assumed  to  be  at  a level  corresponding  to  yield  or  ultimate  stress 
of  the  panel  material  and  is  given  by 
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where 
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yield  or  ultimate  stress 

2m-l-*-c 

2n-l*c 


and  c-0,  0.15  and  O.i  for  S-S,  5-C  and  C-0.  boundarv  conditions,  respec- 
tively. 

The  second  formula  is  associated  with  the  elastic  vave  propagation 
between  integration  points  in  the  short  direction  and  is  given  by 
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For  curved  cylindrical  panels  Che  smallest  At  obtained  from  five 
formulas  is  used  for  the  initial  estimate  for  DEPROP  solutions.  The 
first  two  formulas  are  associated  with  high  frequency  modes  of  cylin- 
drical shells  and  are  given  by 
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In  cases  where  the  length  of  the  panel  in  the  x-direction  is  much 
shorter  than  the  arc  length  in  the  9-direction,  it  was  found  that  the 
flat  panel  formulas  were  more  applicable  than  the  above  three  formulas. 
Thus,  the  fourth  and  fifth  formulas  used  are  based  on  equations  131  and 
132  with  changing  to  and  b to  a0Q  in  equation  131,  and  b to  N to 
M and  C ^ to  in  equation  132. 


The  smallest  time  increment  obtained  from  the  modified  formulas 
used  for  either  flat  or  curved  panels  represents  an  estimated  value  that 
will  generally  give  a stable  solution,  but  does  not  necessarily  repre- 
sent an  optimal  value  for  minimizing  computer  time.  If  a stable  solu- 
tion is  obtained  with  this  estimated  time  increment,  the  time  increment 
can  be  increased  until  two  consecutive  solutions  disagree  significantly 
and  the  penultimate  time  increment  selected  for  future  computer  time 
optimization,  if  desired.  If  the  solution  diverges  using  the  initial 
estimated  time  increment  from  the  formula  procedure,  halving  the  time 
increment  should  easily  result  in  a solution  in  the  stable  range.  The 
DEPROP  program  provides  the  option  of  automatically  using  the  above 
estimated  time  increment  or  having  the  user  select  a value. 


4.2.7  Preload  Static  So^tion 

To  account  for  the  steady-state  airloads  on  the  panel  prior  to 
the  blast  encounter,  the  displacement  components  (u,v,w)  are  determined 
for  a uniform  static  pressure  load  (pQ).  These  displacements  are  used 
as  the  initial  values  to  start  the  dynamic  solution  for  ,:he  transient 
blast  pressure.  The  static  dispxacements  are  obtained  from  all  the 
equations  of  motion  of  the  form  given  in  equation  124  by  a relaxation 


t 


technique  that  permits  a solution  of  these  nonlinear  equations  by 

iteratively  reducing  their  residuals  to  zero.  In  the  case  of  the 

equation  of  motion,  the  residuals  represent  the  accelerations  in  the 

u,v,  and  w directions.  The  initial  trial  values  of  zero  for  the  modal 

displacement  coefficients  u , v , w are  used  to  start  the  relaxation 

mn  mn  mn 

procedure. 

4.2.8  Approximate  Solution  for  Elastic-Plastic  Response  of 
Sandwich  Panels 

The  elastic-plastic  option  of  DEPROP  is  limited  to  handling 
single- layered  isotropic  panels,  so  that  a three-layered  isotropic 
sandwich  panel  is  reduced  within  the  program  to  an  equivalent  single- 
layered panel  based  on  equating  corresponding  extensional  and  bending 
stiffnesses.  A sandwich  panel  with  face  sheets  of  the  same  material 
described  by  a^,  f.  , E and  E^,  can  be  reduced  to  an  equivalent  single- 
layered panel  defined  by  the  following  quantities  in  terms  of  the 
nomenclature  of  figure  40: 


(136a) 


(136b) 


(136c) 


(136d) 
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where 


hg  - equivalent  thickness 

E ■ elastic  modulus  of  face  sheets  (equal  to  aQ/c0) 

Eg  ■ equivalent  elastic  modulus 

E^  - strain  hardening  modulus  of  face  sheets 

g 

E “ equivalent  strain  hardening  modulus 

pg  » equivalent  mass  density 

g 

* equivalent  yield  stress 

4.2.9  Initial  Panel  Imperfections 

There  are  initial  geometric  imperfections  in  actual  panels  on 
an  aircraft  which  should  be  considered  when  a panel  is  being  analyzed  to 
determine  the  threshold  of  permanent  damage.  Generally,  initial  devia- 
tions from  the  cylindrical  shape  are  either  nor.  specified  or  simply  not 
known  during  flight  conditions.  The  displacements  and  stresses  Induced 
in  cylindrical  panels  may  be  sensitive  to  these  imperfections  depending 
on  their  amplitude  and  shape  as  indicated  in  reference  34  for  the  full 
cylindrical  shell . If  initial  radial  imperfection  data  are  available  or 
a nominal  amount  is  specified  as  an  approximation,  such  data  can  be 
included  in  the  analysis  in  the  form  of  modal  imperfection  coefficients, 

A (reference  equation  99). 
mn 

4.2.10  DEPROP  Response  Comparisons  with  Experiment  and  the 
PETROS-3  Structural  Code 

To  evaluate  the  predictive  capabilities  of  DEPROP,  the  calcu- 
lated displacement  and  strain  responses  are  compared  with  available 
experimental  results  (ref.  45)  and  with  results  generated  from  the 
finite  difference  structural  code  PETROS-3  (ref.  46).  The  scope  of 
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these  comparisons  is  not  extensive  enough  to  result  in  a conclusive 
evaluation  of  DEPROP,  but  does  provide  an  initial  evaluation  of  DEPROP 
response  relative  to  available  experimental  data  and  another  structural 
code  based  on  a different  type  of  solution  technique.  For  these  com- 
parison applications  the  DEPROP  routine  of  NOVA  was  temporarily  mated  to 
an  existing  pressure  loading  routine  from  Kaman  AviDyne's  computer 
program  library.  This  special  loading  routine  included  pressure  time 
histories  with  linear  and  exponential  decaying  functions  which  were  used 
as  transient  loadings  for  the  DEPROP  response  comparisons  with  experi- 
ment and  the  PETROS-3  structural  code. 

Some  shock-tube  test  data  were  available  for  panels  from  tests 
performed  at  the  Armament  Laboratory  at  Eglin  Air  Force  Base  and  are 
given  in  reference  45.  The  tests  at  Eglin  essentially  consisted  of 
exposing  clamped  square  panels  to  a blast  from  a disposable  gas-bag 
shock  tube.  The  flat  panels  were  subjected  to  load  levels  that  produced 
large  permanent  deformation  and,  in  some  cases,  rupture  of  the  material 
along  the  clamped  edges.  Primarily,  only  the  permanent  sets  of  the 
deformed  plates  were  measured  and  in  a few  tests  the  deflection  patterns 
were  measured  by  Moire  fringe  photography.  Strain  measurements  were  not 
obtained  for  these  tests.  The  square  flat  panels  are  18  in.  by  18  in. 
and  analytical  and  experimental  comparisons  are  made  for  plate  thick- 
nesses of  0 063  in.  and  0.071  in.  The  panel  material  is  2024-T3  alumi- 
num for  which  a typical  stress-strain  curve  is  shown  in  figure  42  by  the 
solid  curve.  The  bilinear  representation  of  this  stress-strain  curve 
used  in  DEPROP  is  given  by  the  dashed  lines  and  corresponds  to  a yield 
stress  of  50,000  psi  and  a strain  hardening  slope  of  1.24  x 10^  psi. 

The  geometric  and  physical  properties  of  these  plates  are  summarized  in 
table  12.  Pressure  measurements  were  obtained  from  an  instrumented 
rigid  panel  tested  separately  under  essentially  the  same  loading  as  the 
test  panels.  These  pressure  data  were  used  to  determine  an  approximate 
analytical  fit  in  reference  45  using  a uniform  spatial  distribution  with 
the  temporal  decay  given  by 


p(t)=p  (1  — t/r)  e 
m 


(137) 
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TABLE  12 

GEOMETRIC  AND  PHYSICAL  PROPERTIES  OF  PANELS 


Eglin  Teats 

PETROS-3 

Property  of  Panel 

Comparisons 

Comparisons 

Length  (in.) 

18 

20  f 

/ Flat  Panel 

Width  (in.) 

Radius  (in.) 

Subtended  Angle  (Deg) 

18 

20  ) 

100  1 _ . _ . 

, Curved  Panel 

11.48  ! 

Thickness  (in.) 

0.063  and  0.071 

0.1 

Modulus  of  Elasticity  (psi) 

10.5  x 106 

io7 

Poisson's  Ratio 

0.33 

0.3 

Mass  Density  | ^ / 

0.259  x 10~3 

0.259  x IO-3 

Yield  Stress  (psi) 

50,000 

40,000 

Strain  Hardening  Slope  (psi) 

1.24  x 105 

0 and  1.5  x 1G3 

where 


I 


p * 600  psi 

m 

a - 5.7 

t • 1.5  msec 

The  DEPROF  panel  routine  with  the  special  loading  routine  was 
exercised  for  two  test  cases  (panel  thicknesses  of  0.063  and  0.071  in.) 
and  the  predicted  permanent  sets  compared  with  the  test  measurements. 

The  DEPROP  solutions  employed  25  modes  selected  from  a 7x7  array,  a 
15x15  spatial  integration  net  over  the  quarter  panel,  and  five  inte- 
gration points  through  the  thickness.  The  0. 063-in. -thick  panel  under- 
went severe  plastic  deformations  and  was  considered  to  be  near  the 
threshold  of  rupture,  although  it  did  not  fail.  Figure  43  illustrates 
the  analytically  determined  center  displacement  time  history  for  the 
0.063-in.  panel.  This  figure  shows  an  analytical  maximum  displacement 
of  3.43  in.  and  a permanent  set  of  about  3.1  in.  comparing  closely  with 
the  experimental  permanent  set  of  3.0  in.  Figure  44  illustrates  the 
deflection  shape  of  the  panel  at  various  times  during  the  response.  The 
changing  shape  pattern  predicted  by  DEPROP  corresponds  favorably  with 
that  observed  experimentally  through  the  Moire  fringe  technique  and  the 
final  permanent  shape.  Figure  45  shows  the  inner  and  outer  surface 
strain  near  the  center  of  the  panel  where  the  maximum  analytical  strain 
occurred.  The  maximum  analytical  tensile  strain  is  0.133  in. /in.  and 
the  strain  behavior  across  the  thickness  of  the  plate  at  this  position 
is  strongly  membrane.  It  should  be  noted  that  the  maximum  tensile 
strain  probably  occurs  at  the  center  of  the  clamped  edge,  but  the  number 
of  modes  used  analytically  are  not  sufficient  to  predict  the  edge  strain 
accurately.  The  fracture  strain  for  this  material  is  about  0.15  in. /in. 
and,  since  this  panel  did  not  rupture  the  edge  strain,  is  probably  not 
much  greater  than  the  strain  near  the  center  of  the  plate.  For  the 
0.071-in.  panel,  figure  46  shows  the  analytically  determined  center 
displacement  time  history.  Two  panels  of  this  thickness  were  tested 
experimentally  and  the  closeness  of  the  two  permanent  sets  at  2.7  and 
2.8  in.  indicated  good  repeatability  of  the  experiment.  Figure  46 


1 


-xOO  — 


MAX  STRAIN  (IN/IN) 


CENTER  DISPLACEMENT  (IN) 


indicates  that  DEPROP  predicts  a permanent  set  of  about  2.65  in.  which 
compares  well  with  the  experimental  values.  The  differences  in  the 
analytical  and  experimental  permanent  set  values  for  -he  two  test 
plates  considered  are  less  than  4%. 

Comparisons  of  displacement  and  strain  responses  from  DEPROP 
with  PETROS-3  are  made  for  various  panels  subjected  to  a simple  trans- 
ient pressure  loading.  PETROS-3  (ref.  44)  is  a dynamic  response  struc- 
tural code  for  the  calculation  of  large  elastic-plastic  deformations  of 
plates  and  shells  based  on  the  finite-difference  method  of  solution. 
While  DEPROP  is  based  on  a modal  type  solution,  there  are  several  simi- 
larities between  the  numerical  methods  used  in  DEPROP  and  PETROS-3. 

Both  codes  use  the  same  temporal  numerical  method  (central-difference) 
and  the  same  Gaussian  integration  technique  through  the  thickness  of  the 
panel  for  elastic-plastic  solutions.  It  should  be  noted  that  neither 
DEPROP  nor  PETROS-3  can  be  considered  as  an  absolute  standard  in  these 
comparisons.  Which  solution  is  more  accurate  in  computing  displacement 
and  strain  quantities  within  a panel  can  only  be  determined  by  thorough 
correlation  with  well  controlled  experiments.  Solutions  obtained  using 
DEPROP  and  PETROS-3  were  based  on  the  respective  mode  or  mesh  limits 
presently  dimensioned  in  each  program  for  reasonable  computer  running 
times  and  core  size  on  the  CDC  6600  computer.  Thus,  the  solutions  are 
not  necessarily  optimal  relative  to  complete  convergence  of  strain 
quantities  throughout  the  panels.  For  example,  neither  the  number  of 
modes  used  in  the  DEPROP  solutions  nor  the  finite-difference  mesh  size 
used  in  the  PETROS-3  solutions  are  sufficient  to  determine  accurately 
the  strains  at  the  clamped  edges  of  a panel  undergoing  large  plastic 
deformations  where  large  strain  gradients  exist  very  near  the  edge.  For 
this  reason,  only  comparisons  are  made  with  center  strains  of  the  panels. 
For  all  the  panel  solutions  using  DEPROP  and  PETROS-3,  an  adequate  time 
increment  was  selected  so  that  numerical  convergence  was  achieved  on  a 
temporal  basis  and  the  same  number  of  integration  points  through  the 
thickness  was  used  in  both  code  solutions.  The  panel  problems  selected 
were  all  symmetrical  so  that  only  one  quarter  of  the  panel  was  con- 
sidered in  all  solutions.  For  DEPROP  solutions,  35  modes  were  selected 
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from  a 7x7  array  and  a 15x15  spatial  integration  net  was  used  on  the 
quarter  panel.  For  PETROS-3  solutions,  a uniform  finite-difference  mesh 
of  16x16  on  the  quarter  panel  was  modified  by  halving  the  mesh  size  near 
the  center  and  edges  of  the  panel  to  produce  a nonuniform  20x20  mesh. 
Comparisons  of  center  displacement  and  strain  responses  from  DEPROP  and 
PETROS-3  were  made  for  simply  supported  and  clamped  20  in.  x 20  in. 
x 0.1  in.  square  panels  with  purely  elastic,  elastic-perf ectly  plastic 
and  elastic-strain  hardening  materials.  A clamped  rectangular  panel  and 
a clamped  cylindrical  panel,  both  with  material  strain  hardening,  were 
also  used  for  the  response  comparisons.  The  geometric  and  physical 
properties  of  the  various  panels  are  summarized  in  table  12.  Through 
the  special  loading  routine,  these  panels  were  subjected  to  a uniform 
pressure  loading  with  a triangular  pulse  shape  that  instantaneously  rose 
to  the  peak  pressure  at  zero  time  and  decayed  to  zero  pressure  at  2 
msec . 

The  DEPROP  and  PETROS-3  comparisons  of  center  displacement  and 
strain  responses  for  various  material  conditions  are  shown  in  figures  47 
through  52  for  the  simply  supported  square  panels  at  a peak  pressure  of 
200  psi;  in  figures  53  through  59  for  the  clamped  square  panels  at  peak 
pressures  of  250  psi  (elastic-plastic  solutions)  and  100  psi  (purely 
elastic  solution);  in  figures  60  and  61  for  the  clamped  rectangular 
panel  (aspect  ratio  of  1.5)  at  a peak  pressure  of  250  psi;  and  in 
figures  62  through  64  for  the  clamped  cylindrical  panel  at  a peak 
pressure  of  250  psi.  From  these  comparisons  of  DEPROP  and  PETROS-3 
responses,  the  following  observations  are  made: 

1.  The  comparisons  for  purely  elastic  large  displace- 
ment response  of  a square  clamped  panel  (figures  53 
and  54)  show  satisfactory  agreement  between  DEPROP 
and  PETROS-3  for  both  center  displacements  and 
strains . 

2.  For  large  elastic-plastic  deformations  of  these 
panels,  the  comparisons  indicate  that  the  best 
agreement  (magnitude  differences  less  than  5J») 
occurs  for  all  center  displacement  responses 
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Comparison  of  DEPROP  and  PETROS- 3 Displacement  Resoonse 
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Figure  3J.  ZocDanscn  j:  DIPRC?  ir.d  PET! 
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Figure  62.  Comparison  of  DEDROP  and  PETROS-3  Displacement  Response  j 

for  a Clamped  Curved  Panel  3 
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Figure  63.  Comparison  of  DEPROP  and  PETROS-3  Axial  Strain  Response 
for  a Clamped  Curved  Panel. 
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and  strain  response  in  the  simply  supported 
and  perfectly  plastic  plate  (figures  47-49); 
but  as  the  clamped  edges,  strain  hardening  or 
curved  geometry  are  introduced,  the  differences 
in  strain  magnitude  increase,  generally  varying 
betweeu  about  5 and  25%. 

3.  The  general  behavior  of  the  DEPROP  and  PETROS-3 
response  solutions  is  the  same  throughout  the 
panels,  the  primary  differences  are  characterized 
by  larger  peak  magnitudes  and  times  of  peak 
response  being  predicted  by  PETROS-3  for  both 
displacement  and  strain  responses.  The  largest 
differences  seem  to  be  associated  with  the  clamped 
boundary  conditions.  Near  the  end  of  this  inves- 
tigation, PETROS-3.5  (ref.  46)  became  available. 
PETROS-3.5  is  a revised  version  of  PETROS-3  that 
has  improved  the  finite-difference  representation 
near  the  boundaries  and  the  numerical  scheme  used 
in  the  Incremental  plasticity  model.  The  elastic- 
plastic  solution  for  the  square  clamped  panel  was 
obtained  using  PETROS-3.5  and  compared  with  the 
PETROS-3  results.  It  was  found  that  differences 
in  the  central  displacement  and  strain  responses 
were  less  than  2%  between  the  PETROS-3  and  3.5 
solutions  for  the  selected  panel. 


4.  At  early  times  during  the  response,  the  DEPROP 

center  strain  response  exhibits  some  oscillations 
not  produced  in  the  PETROS-3  strain  response.  In 
the  PETROS-3  solution  at  early  times,  the  central 
position  of  the  plate  remains  absolutely  flat, 
exhibiting  just  membrane  strains  in  this  region  of 
the  plates.  However,  based  on  the  modal-type 
solution  in  DEPROP,  this  central  portion  of  the 
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place  is  only  nearly  flat  and  exhibits  bending 
strains  as  well  as  membrane  strains  in  this  region 
which  account  for  the  oscillations  in  the  strain 
plots  at  early  times. 

5.  It  was  found  that  the  convergence  of  the  strain 
response  in  the  PETROS-3  calculation  is  much  more 
sensitive  to  increases  in  time  increment  and  mesh 
size  than  the  DEPROP  calculation  is  to  Increases  in 
time  Increment  and  decreases  in  the  number  of  modes 
used. 

6.  Because  of  the  larger  strain  gradient  near  the 
clamped  edges  of  the  severely  deformed  panels,  it  is 
felt  that  neither  DEPROP  nor  PETROS-3  accurately 
predicts  the  strains  at  a clamped  edge  with  the 
modes  and  finite-difference  mesh,  respectively,  used 
herein.  To  assess  the  degree  of  accuracy  of  analy- 
tical predictions,  experimental  strain  data  are 
needed  near  the  clamped  edges  or  selected  panels 
loaded  with  a well  defined  pressure  time  history. 

7.  Although  neither  code's  solution  parameters,  such  as 
modes  and  mesh  size,  should  be  considered  as  optimal 
for  the  computation  of  converged  center  strains,  it 
was  found  that  PETROS-3  required  a smaller  time 
increment  than  DEPROP  for  a nearly  convergent 
center  strain  solution.  For  these  panels  using  the 
aforementioned  solution  parameters,  DEPROP  solutions 
used  about  half  the  central  processor  computer  time 
used  by  PETROS-3  solutions. 
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SECTION  V 
DAMAGE  CRITERIA 


The  computer  program  contains  the  option  of  selecting  two  levels  of 
damage  to  the  structural  elements  being  analyzed.  One  level  corre- 
sponds to  no  permanent  damage  and  the  other  level  corresponds  to  catas- 
trophic damage.  However,  the  analyst  must  Interpret  the  effect  that  the 
damaged  panel  or  structural  element  has  on  the  performance  or  mission 
of  the  aircraft.  The  analyst  also  has  the  option  of  specifying  the  prob- 
ability that  the  damage  level  selected  will  be  exceeded.  For  example, 
the  analyst  might  specify  the  no-damage  level  and  a probability  of  5 
percent  that  the  no-damage  level  will  be  exceeded;  that  is,  that  damage 
will  occur.  The  method  of  establishing  the  probability  that  a specific 
level  of  damage  will  be  exceeded  will  now  be  described. 

Designate  the  structural  response  parameter  (stress  or  strain) 

by  R and  the  value  of  R at  which  the  specified  damage  level  occurs  by 

Rj.  The  preblast  value  of  the  parameter  will  be  designated  by  Rq  and 

the  maximum  value  calculated  by  the  program  with  respect  to  time  by  R^. 

The  value  of  R will  be  defined  such  that  the  probability  of  exceeding 
P 

R,  will  be  the  specified  value,  m.  This  is  accomplished  by  estimating 
a 

the  accuracy  of  the  prediction  and  then  assuming  that  the  probability 
density  distribution  of  the  response  is  normal.  The  assumption  of  a 
normal  distribution  is  based  upon  the  central-limit  theorem,  which 
states  that  the  sum  of  independent  variates  from  the  same  or  different 
distributions  is  normally  distributed  in  the  limit  and  that  this  .limit 
is  approached  very  rapidly  (see  ref.  47).  The  large  number  of  variables 
which  influence  the  response  of  the  structure  justifies  the  assumption 
of  a normal  distribution. 
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Assuming  that  the  analysis  Is  unbiased,  R is  the  mean  of  the 
normal  distribution  illustrated  in  figure  65.  The  saaded  area  repre- 
sents the  probability  that  the  specified  damage  level  will  be  exceeded 
by  m percent.  It  is  assumed  that  the  Inaccuracy  in  Rp  can  be  defined  by 
a factor  X,  such  that  three  standard  deviations  in  response  will  be 
equal  to  XRp;  that  is, 


3o  - XR 


(138) 


The  choice  of  three  standard  deviations  is  arbitrary.  The  problem  of 
defining  X so  that  it  corresponds  to  the  number  of  standard  deviations 
choset  will  be  treated  later. 

The  number  of  standard  deviations,  n,  associated  with  the  prob- 
ability, m,  can  be  determined  from  normal  probability  distribution 
tables.  For  example,  if  m is  5 percent,  the  associated  value  of  n is 
1.645.  Assembling  all  of  the  above  assumptions  and  observations,  the 

value  of  R sought,  which  we  will  denote  by  R , can  be  written  as 
P P 


R ■ R , - no 
P d 


(139) 


or,  introducing  equation  138  with  the  response,  Rp,  replaced  by  the 
desired  response,  Rp, 


R - R , - n | R 
p d 3 p 


(140) 


Solving  equation  140  for  Rp, 


R m , , . 
P 1 + <; 


(141) 
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PROBABILITY  DENSITY  FUNCTION 


Rp  Rp 

RESPONSE  PARAMETER,  R 

*Provided  that  R = r 

P P 


Figure  65  - Normal  Probability  Density  Distribution  of 
Response  Parameter,  R 


i 
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where 


5 = y (142) 

The  remaining  problem  In  Che  implementation  of  the  above  equations  is 

the  determination  of  the  structural  response,  R , , associated  with  the 

d 

level  of  damage  specified,  and  the  determination  of  the  inaccuracy  fac- 
tor, X. 


5.1  DETERMINATION  OF  DAMAGE  PARAMETER,  Rj 

d 

The  program  contains  the  option  of  selecting  two  damage  levels, 
no-damage  and  catastrophic  damage.  The  no-damage  response  level  is 
Interpreted  as  that  level  of  structural  response  which  represents  the 
threshold  of  damage.  The  severe  or  catastrophic  level  of  damage  is 
interpreted  as  that  level  of  response  at  which  the  material  ruptures  or 
fractures. 

When  the  damage  level  and  the  percent  probability  that  the  level 
will  be  exceeded  are  specified,  the  program  can  be  used  to  determine 
the  range  at  which  the  desired  response  will  occur  with  the  desired 
probability.  The  process  is  iterative;  a trial  range  is  specified,  and 
the  response  is  determined  and  compared  to  a criterion.  The  range  is 
adjusted  within  the  program  and  the  process  is  repeated  until  the  ratio 
of  maximum  stress  (Rp)  to  the  critical  stress  (Rp)  converges  to  unity. 

For  no-damage  response,  the  program  uses  the  largest  value  f the 

ratio  of  maximum  stress  (R  ) to  the  critical  stress  (R  ) . In  order  to 

P P 

permit  an  orderly  iteration  process,  values  of  the  ratio  (R  /R  ) greater 

than  unity  must  be  allowed.  In  addition,  the  critical  stress  may  be 
greater  than  the  yield  stress,  depending  upon  the  selected  value  of 
probability,  m.  Actually,  the  material  yields  and,  presumably,  the  struc- 
ture follows  a different  branch  of  the  stress-strain  curve;  hence,  an 
artificial  means  must  be  incorporated  within  the  program  to  permit 
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hi'/’.ier  stresses  than  the  yield  stress  along  the  elastic  curve.  This  is 
accomplished  by  disregarding  yielding  and  using  the  extended  elastic 
stress-strain  curve  for  the  no-damage  condition.  This  artifice  permits 
a continuous  variation  of  the  response  with  range,  which  is  necessary 
in  order  to  achieve  a systematic  iteration  process. 

The  process  is  similar  for  catastrophic  damage.  For  plastic  materials 
or  for  buckling  analyses,  the  process  is  identical,  except  that  larger 
stresses  are  considered  critical.  Otherwise,  the  elastic-plastic  stress- 
strain  curve  must  be  made  to  accommodate  strains  exceeding  the  fracture 
level.  In  DEPROV,  the  stress-strain  curve  is  extended  beyond  ultimate 
strain  with  the  same  strain  hardening  slope,  while  in  DEPROB  a perfectly 
plastic  extension  of  the  stress-strain  curve  is  assumed.  Both  of  these 
artificial  extensions  permit  the  calculation  of  response  strains  beyond 
the  ultimate  strain  and  therefore  allows  a smooth  iteration  process 
for  all  selected  values  of  probability,  m. 

The  type  of  structure  being  analyzed  determines  the  level  of  stress 
or  deformation  associated  with  the  threshold  of  damage  or  catastrophic 
damage.  Both  levels  of  damage  are  discussed  for  single-layered  and  multi- 
layered honeycomb  panels  and  for  stiffeners,  frames,  radomes,  and  ribs. 

5.1.1  Single-Layered  Panels 

The  DEPROP  program  is  used  to  determine  the  stresses  and  strains 
induced  in  single-layered  panels  by  steady-state  and  transient  pressure 
loads.  These  stress  and  strain  quantities  are  related  to  conditions  in 
the  material  that  can  produce  permanent  deformations  or  rupture  of  a 
panel.  Yielding  of  the  material  is  taken  as  the  limiting  condition 
for  threshold  of  permanent  damage  for  metal  panels.  Since  the  panel 
deformations  are  biaxial,  the  yield  stress  is  compared  with  the  equiv- 
alent stress  a associated  with  the  Mises-Hencky  yield  criterion 
(equation  109).  This  yield  criterion  states  that  plastic  flow  will 
occur  when  the  equivalent  stress  a reaches  a value  equal  to  the  uniaxial 
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yield  stress.  For  plastic  or  reinforced  plastic  materials  which  are 
assumed  to  be  loaded  elastically  to  an  ultimate  stress  or  strain  and 
therefore  attain  a fracture  condition  without  actually  yielding,  70 
percent  of  the  tensile  ultimate  stress  and  100  percent  of  the  compressive 
ultimate  stress  are  taken  as  the  critical  stresses  which  define  the 
limiting  condition  for  threshold  of  permanent  damage.  It  should  be 
noted  that  curved  panels  could  undergo  elastic  buckling  prior  to  material 
yielding,  but  this  phenomenon  is  not  considered  critical  since  after  the 
load  Is  removed  the  panel  returns  to  its  original  undeformed  condition. 

For  catastrophic  damage  of  single-layered  metal  panels,  large 
inelastic  deformations  are  produced  during  the  response  that  lead  to 
rupture  of  the  material.  An  approximate  rupture  criterion  is  established 
which  compares  the  uniaxial  rupture  fracture  strain  of  the  material 
with  the  effective  strain  e (equation  112)  using  checks  to  insure  that 
a tensile  strain  condition  is  present.  For  catastrophic  damage  of  plas- 
tic  and  reinforced  plastic  panels,  the  criterion  is  based  on  the  ulti- 
mate tensile  stress  of  the  material  and  is  compared  with  the  principal 
tensile  stresses  in  the  panel.  For  brittle  materials,  it  is  assumed 
that  fracture  occurs  at  the  ultimate  tensile  stress. 


5.1.2  Honeycomb  Panels 

The  multilayer  panel  option  of  the  DEPROP  program  is  used  to 
determine  elastic  stresses  and  strains  induced  in  honeycomb  panels  due 
to  pressure  loads.  The  honeycomb  panel  is  a three-layered  panel  with 
either  isotropic  or  orthotropic  material  properties. 

The  yield  stress  in  the  fr.ee  sueet  is  used  as  one  criterion 
or  limiting  condition  for  metal  face  sheets  when  establishing  the 
threshold  of  permanent  damage  for  honeycomb  panels.  For  reinforced 
plastic  face  sheets,  70  percent  of  tensile  ultimate  stress  and  100  per- 
cent of  compressive  ultimate  stress  are  used  for  the  threshold  of 
permanent  damage.  Other  criteria  for  honeycomb  panels  are  often  required 
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since  the  panels  can  experience  local  instabilities  such  as  intracell 
buckling  or  face  sheet  wrinkling.  The  intracell  buckling  occurs  mainly 
for  honeycomb  panels  with  thin  face  sheets,  and  face  sheet  wrinkling 
occurs  mainly  for  low-density  cores  and  weak  bonding  systems.  If  these 
instabilities  are  severe  enough  to  result  in  permanent  deformations  in 
the  panel,  they  are  acceptable  as  a Uniting  condition  for  threshold 
of  permanent  damage.  Formulas  are  given  in  reference  48  which  relate 
intracell  buckling  or  wrinkling  to  the  geometric  and  material  proper- 
ties cf  the  honeycomb  panel  as  follows.  For  intracell  buckling, 


a - 3E, 
cr  f 


(143) 


a - critical  stress  in  the  face  sheet,  psi 
cr  r 

f ■ face  sheet  thickness,  inches 

d » core  cell  size,  inches 

Ef  ■ elastic  modulus  of  the  face  sheet,  psi 


For  wrinkling. 


a - 0.5  (G  E EJ1/3 
cr  c c f 


(144) 


* critical  stress  in  the  face  sheet,  psi 

■ core  shear  modulus,  psi 

■ core  modulus  of  elasticity  parallel  to  the  core  depth,  psi 


modulus  of  elasticity  of  the  face  sheet,  psi 
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As  indicated  earlier,  the  elastic-plastic  option  of  DEPROP  is 
limited  to  single-layered  isotropic  panels.  Consequently,  in  deter- 
mining catastrophic  damage  of  honeycomb  metal  panels,  an  "equivalent" 
single-layered  panel  is  developed  in  subsection  4.2.8  to  represent  the 
deformation  response  of  the  original  honeycomb  in  the  elastic  and 
inelastic  ranges.  Catastrophic  damage  for  a honeycomb  metal  panel  is 
based  on  the  rupture  strain  of  the  face  sheet  material. 

5.1.3  Stiffeners  and  Frames 

The  DEPROB  routine  can  be  used  to  calculate  the  response  of 
stiffeners  and  frames  combined  with  local  effective  skin  due  to  a 
pressure  loading.  The  dynamic  response  includes  deflections  and  accel- 
erations of  the  structural  element  plur  stress  throughout  the  cross 
section,  particularly  stress  in  the  outer  fibers  of  the  cross  section. 
The  stresses  of  the  outer  fibers  apply  to  either  the  flange  plus  the 
effective  skin  attached  to  the  flange  or  to  the  outstanding  leg  of 
the  element.  In  general,  either  can  be  loaded  in  tension  or  com- 
pression. To  determine  the  threshold  of  permanent  damage,  the  stresses 
in  the  outer  fibers  are  compared  to  the  yield  stress  for  tensile  or 
compressive  loads.  The  magnitude  of  the  moment  at  the  fixed  ends  of  a 
uniform  beam  carrying  a uniform  static  load  is  twice  that  of  the  moment 
at  t.ie  center  of  the  beam.  Consequently,  the  threshold  of  permanent 
damage  will  be  achieved  primarily  by  buckling  of  the  outstanding  leg  at 
che  eud  of  the  beam  for  this  case,  rather  than  yielding  in  tension  or 
compression  at  any  point  along  the  beam.  Crippling  stress  formulas  are 
available  for  local  buckling  of  outstanding  legs  of  different  shaped 
stringer  or  frame  elements.  The  crippling  stress  formula  used  in  the 
program  is  identical  to  that  used  in  reference  49,  page  Cl. 2-36: 
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where 


f - crippling  stress  for  the  outstanding  leg,  psi 

K - buckling  coefficient  equal  to  1.25  if  the  outstanding  leg 

has  only  one  corner,  and  equal  to  3.62  if  the  outstanding  leg 
has  two  corners 

E - modulus  of  elasticity,  psi 

t^  - thickness  of  the  outstanding  leg,  inches 

w - width  of  the  outstanding  leg,  inches 

For  nonuniform  beam  elements,  t^  and  w at  the  center  of  the  beam  are 
used  as  representative  values. 

Local  buckling  or  crippling  is  a minor  type  of  failure  com- 
pared with  reaching  ultimate  strain  or  rupture  of  some  portion  of  the 
cross  section.  Consequently,  catastrophic  damage  for  stiffeners  and 
frames  is  based  on  tensile  rupture  strain  in  the  outer  fibers. 

And  for  clamped  beam  elements  subjected  to  catastrophic 
damage,  DEPROB  computes  a very  localized  ideal  edge  strain  which  can  far 
exceed  usual  handbook  rupture  strain  levels  prior  to  actual  rupture. 
Therefore,  a strain  equal  to  one-third  that  computed  at  the  edge  is 
compared  with  rupture  strain  levels  for  catastrophic  damage. 

5.1.4  Radomes  and  Other  Shells 

Radomes  on  various  aircraft  have  different  shapes.  Some 
radomes  are  best  analyzed  by  the  DEPROP  program  where  a curved  panel 
representation  is  reasonable.  Other  radomes,  such  as  the  nose  and  tail 
radomes  of  the  Bl,  are  conical  or  near-cylindrical  shells.  For  these 
shapes,  a two-dimensional  ring  representation  is  reasonable  and  the 
DEPROB  program  should  be  used.  Since  reinforced  plastic  material  which 
fractures  at  ultimate  strain  is  used  in  radomes,  70  percent  of  tensile 
ultimate  stress  and  100  percent  compressive  ultimate  stress  are  used  fcr 
the  threshold  of  permanent  damage,  and  ultimate  tensile  strain  is  used 
for  catastrophic  damage. 
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DEPROB  can  also  be  used  to  analyze  a metal  shell  In  an  approx- 
imate way  when  modeled  as  a free  ring.  In  this  case  the  threshold  of 
damage  criterion  Is  based  on  tensile  and  compressive  yield  stresses, 
compared  to  stresses  In  the  outer  fibers.  Catastrophic  damage  Is  based 
on  tensile  strain  In  the  outer  fibers  exceeding  rupture  strain. 

5.1.5  Rib  Webs  and  Stiffeners 

Failure  In  vertical  metal  rib  elements  is  assumed  to  occur  In 
a buckling  mode.  Modeled  as  an  axially  loaded  beam  with  fixed  ends 
(except  that  motion  is  permitted  in  the  axial  direction),  threshold  of 
permanent  damage  and  catastrophic  damage  are  both  related  to  the  occur- 
rence of  column  buckling.  It  should  be  noted  that  as  a material  begins 
to  yield,  it  quickly  loses  its  resistance  to  buckling.  Although  the  two 
events  are  clearly  separate  phenomena,  experience  indicates  that  if  for 
Increasing  levels  of  load  the  beam  has  not  buckled  by  the  time  it  begins 
to  yield,  it  will  with  very  little  additional  load. 

Catastrophic  damage,  then,  is  defined  as  the  point  at  which 
the  tensile  strain  in  the  outer  fibers  exceeds  the  tensile  ultimate 
strain,  which  will  usually  closely  follow  the  attainment  of  tensile 
yield  strains.  Threshold  of  permanent  damage  is  defined  as  a maximum 
tensile  or  compressive  strain  equal  to  70Z  of  yield  strain. 

For  situations  in  which  the  rib  collapses,  the  numerical  pro- 
cedures employed  may  not  be  able  to  keep  up  with  the  process.  In  such 
cases,  the  program  assigns  a negative  number  to  the  maximum  response  and 
proceeds  to  select  a larger  range. 

5.1.6  Damage  Criteria  Summary 

In  summary,  damage  criteria  for  panels,  stiffeners,  frames, 
radomes,  and  ribs  are  given  in  table  13  for  threshold  of  permanent 
damage,  T?D,  and  catastrophic  damage,  CD.  In  general,  all  spatial 
locations  on  the  structural  elements  are  checked  in  the  timewise  solu- 
tion in  order  to  accumulate  the  maximum  response  parameters.  While 
DEPROB  checks  the  parameters  every  time  step,  DEPROP  checks  only  every 
tenth  step-a  procedure  which  conserves  computer  time  in  the  slower 
running  DEPROP  without  significant  loss  of  accuracy. 


i 
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TABLE  13 
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CRITERIA  FOR  THRESHOLD  OF  PERMANENT  DAMAGE  (TPD) 
AND  CATASTROPHIC  DAMAGE  (CD) 


Program  - DEPROP 
Structural  Element  - Panels 

Spatial  Locations  Considered  - All  spatial  Integration  points 

Types  of  Panels: 

1.  Single-layered  metal 

TPD  - Compare  tensile  or  compressive  yield  stress  with 
Mises-Hencky  equivalent  stress  at  inner  and  outer 
surfaces . 

CD  - Compare  rupture  strain  with  effective  tensile  strain 
at  inner  and  outer  surfaces. 


2.  Single-layered  plastics 

TPD  - Compare  70  percent  of  tensile  and  100  percent  of 

compressive  ultimate  stresses  with  principal  tensile 
and  compressive  stresses,  respectively,  at  the  inner 
and  outer  surfaces. 

CD  - Compare  tensile  ultimate  stress  with  principal  ten- 
sile stresses  at  inner  and  outer  surfaces. 


3.  Honeycomb  metal 

TPD  - Compare  tensile  yield  stress  with  Mises-Hencky  equiva- 
lent stresses  at  centers  of  face  sheets;  compare  the 
compressive  yield  stress  with  Mises-Hencky  equivalent  ! 
stress  at  centers  of  face  sheets;  and  compare  the  low- 
est of  intracell  buckling  stress  and  face  sheet  wrin- 
kling stress  with  maximum  principal  compressive 
stresses  at  centers  of  face  sheets. 
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TABLE  13  (Continued) 
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CD  - Compare  face  sheet  rupture  strain  with  effective 
tensile  strains  at  points  on  equivalent  single* 
layered  cross  section  which  correspond  to  the  center 
of  the  face  sheets  of  original  cross  section. 

4.  Honeycomb  reinforced  plastics 

TPD  - Compare  70  percent  of  tensile  ultimate  stress  with 

principal  tensile  stresses  at  centers  of  face  sheets; 
and  compare  the  lowest  of  100  percent  compressive 
ultimate  stress,  intracell  buckling  stress,  and  face 
sheet  wrinkling  stress  with  maximum  compressive 
stresses  at  centers  of  face  sheets. 

CD  - Compare  ultimate  tensile  stress  with  principal 
tensile  stresses  at  centers  of  face  sheets. 


5.  Multilayered  plastics 

TPD  - Compare  70  percent  of  tensile  and  100  percent  of 

compressive.  ultimate  stresses  with  principal  tensile 
and  compressive  stresses,  respectively,  at  inner  and 
outer  surfaces  of  each  layer. 

CD  - Compare  ultimate  tensile  stress  with  principal 

teusile  stresses  at  inner  and  outer  surfaces  of  each  | 
layer. 

i 

Program  - DEPROB 

I 

Structural  Element  - Stiffeners  and  Frames  (beam  elements) 

Spatial  Locations  Considered  - All  bars  of  beam 

TPD  - More  critical  of  (a)  comparison  of  largest  tensile 
and  compressive  stresses  on  inner  and  outer  fibers 
with  tensile  and  compressive  yield  stresses,  respec-  , 
tively;  (b)  comparison  of  largest  compressive  stress 
on  outstanding  leg  with  the  crippling  stress. 

CD  - Compare  rupture  strain  with  largest  tensile  strains  j 
on  inner  and  outer  fibers.  j 


l 

i 

l 
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TABLE  13  (Concluded) 


! 

Structural  Element  - 

Shells  (Free  Ring) 

Spatial  Locations  Considered  - All  Bars  of  Ring 

• 

l 

[ ( 

Types  of  shells 

: 

; 

; 

1.  Reinforced  Plastic  Radomes 

[ 

t 

\ 

TPD 

- Compare  70  percent  of  tensile  and  100  percent 

of  compressive  ultimate  stresses  with  tensile 
and  compressive  stresses,  respectively,  at 
inner  and  outer  fibers. 

[ 

CD 

- Compare  ultimate  tensile  strain  with  largest 

tensile  strains  on  inner  and  outer  fibers. 

\ 

. 

2.  Metal 

ii 

i 

i 

r 

; 

TPD 

Compare  tensile  and  compressive  yield  stresses 
with  tensile  and  compressive  stresses,  respec- 
tively, in  the  inner  and  outer  fibers  of  each 
layer . 

: 

> 

CD 

Compare  rupture  strain  with  largest  tensile 
strains  on  inner  and  outer  fibers. 

‘ 

Structural  Element  - 

Ribs  (End  Loaded,  Clamped  Beam  Element) 

Spatial  Locations  Considfved  - All  Bars  of  Beam 

: 

. 

' 

l.r’D 

Compare  70  percent  of  tensile  and  compressive 
yield  stress  with  largest  tensile  and  com- 
pressive stresses,  respectively,  in  inner  and 
outer  fibers. 

1 

CD 

Compare  rupture  strain  with  largest  tensile 
strains  on  inner  and  outer  fibers. 

5.2  DETERMINATION  OF  INACCURACY  FACTOR,  X 

In  order  to  establish  the  probability  of  threshold  or  catastrophic 
damage  occurring,  the  inaccuracy  factor,  X,  in  equation  (142)  must  be 
defined.  Note  tiiat  the  determination  of  X is  highly  subjective  and  is 
not  subject  to  verification  by  analysis.  Only  by  comparing  a large 
group  of  experimentally  produced  responses  to  pressure  loads  from 
nuclear  blasts  under  varying  conditions  with  corresponding  calculated 
responses  could  X be  determined  objectively. 

Assuming  that  the  environmental  conditions  are  known  exactly,  the 
largest  Inaccuracies  are  assumed  to  exist  in  the  modeling  of  the  struc- 
ture as  discrete  elements.  Values  of  fractional  inaccuracy,  X,  have 
been  selected  for  the  various  conditions  considered  and  related  to 
accuracy  factors  which  are  given  in  table  14.  These  values  are  very 
subjective  and  represent  the  best  estimates  the  authors  can  make  based 
on  their  experience  in  the  present  and  related  problems.  If  the  accuracy 
factor  is  2,  the  actual  response  is  considered  to  be  within  a factor 
of  2 of  the  predicted  response.  This  results  in  different  values  of  X 
(and  hence  standard  deviations)  for  responses  greater  than  Rp(X+)  and 
less  than  R^  (X  ) (figure  65).  These  values  cf  X are  related  to  ACC  by 


ACC  - 1 


(146) 


1 - 


ACC 


(147) 


In  terms  of  probability,  if  the  probability,  m,  < 0.5,  X is  used  in  the 
program  to  determine  £;  if  m > 0.5,  X is  used.  The  accuracy  factor  is 
considered  to  be  more  meaningful  than  a single  value  of  X.  For  example, 
if  X were  taken  as  one,  there  would  be  an  implication  that  the  likeli- 
hood of  zero  response  would  correspond  to  the  likelihood  of  twice  the 
predicted  response.  Such  an  implicatijn  can  be  rejected  on  purely 
intuitive  grounds;  hence,  the  accuracy  factor,  which  avoids  this  diffi- 
culty, seems  to  be  more  logical. 
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